Exploring metabolism-growth relationships within and across individual fish.
E. Hoots and P. Biro
2025-08-26
READ ME
This R code is used to take raw respirometry slopes measured on Galaxias maculatus at Deakin’s Marine Research Centre in Queenscliff, Victoria, Australia and calculate their standard, routine, and maximum metabolic rate. These metabolic measurements are made repeatedly on known individuals held under one of two temperature treatments (18°C or 23°C), and paired with repeated mass and length measurements over five month of data collection (May - September, 2023). After reading in all data, a multivariate mixed effects modeling approach is undertaken in order to assess covariance between growth and metabolism traits among- and within-individuals.
CONTACT
Elizabeth C. Hoots
Email: beth.hoots@research.deakin.edu.au
GitHub: https://github.com/Bhoots99
LOAD REQUIRED PACKAGES
# installs and loads packages, need to install "pacman" first
pacman::p_load("brms",
"broom",
"broom.mixed",
"car",
"coda",
"dplyr",
"ggeffects",
"ggiraphExtra",
"ggplot2",
"ggpubr",
"hms",
"lme4",
"lmerTest",
"lubridate",
"nlme",
"posterior",
"readxl",
"tidyverse"
)Working directories
wd <- getwd() # getwd tells us what the current wd is, we are using this to drop it in a variable called wd
raw_data_wd <- paste0(wd, "./00_raw_data") # creates a variable with the name of the wd where raw data are stored
MRcalcs_wd <- paste0(raw_data_wd, "./MRcalcs")
MR_slopes_wd <- paste0(raw_data_wd, "./MR_Slopes")
MMR_wd <- paste0(raw_data_wd, "./MMR")
metadata_wd <- paste0(wd, "./01_metadata") # creates a variable with the name of the wd where metadata are stored
figures_wd <- paste0(wd, "./03_figures") # creates a variable with the name of the wd where figures are to be storedRead in data, calculate metabolic rates, merge into master data table
#calcSMR function adapted from Chabot, 2016 (https://doi.org/10.1111/jfb.12845)
#note that we have not used the mlnd function as the mean of the lowest 10th percentile was better suited to our data
calcSMR <- function(Y) {
u <- sort(Y)
tenpc <- round(0.1 * length(u))
SD10pc <- sd(u[1:tenpc])
low10pc = mean(u[(which((u > (mean(u[1:tenpc])-SD10pc)))):(tenpc+which((u > (mean(u[1:tenpc])-SD10pc -u[1]))))])
return(list(low10pc = low10pc))
}
# Define a function to extract date strings from file names
extract_date <- function(files) {
gsub(".*?([0-9]{8}).*", "\\1", basename(files))
}
# Specify the directories where your files are located
csv_folder <- MR_slopes_wd
xlsx_folder <- MMR_wd
# Get the list of files in each folder
csv_files <- list.files(csv_folder, pattern = "_MR_slopes.csv", full.names = TRUE)
xlsx_files <- list.files(xlsx_folder, pattern = "_MMR.xlsx", full.names = TRUE)
# Extract date strings from file names
csv_dates <- extract_date(csv_files)
xlsx_dates <- extract_date(xlsx_files)
matching_files <- Map(function(date) {
csv_file <- csv_files[grep(date, csv_dates)]
xlsx_file <- xlsx_files[grep(date, xlsx_dates)]
return(list(csv = csv_file, xlsx = xlsx_file))
}, unique(c(csv_dates, xlsx_dates)))
# Loop through matching files and perform your data processing
for (i in seq_along(matching_files)) {
date <- unique(csv_dates)[i]
current_files <- matching_files[[i]] # Renamed 'files' to 'current_files'
# Process CSV data
tb_respirometry <- read_csv(current_files$csv) %>%
rename(
chamber_ch = Chamber.No,
ID_fish = Ind,
mass = Mass,
length = Length,
volume_ch = Ch.Volume,
DOunit = DO.unit,
dateTime = Date.Time,
Temp_class = Temp.class,
slope_wBR = Slope.with.BR,
BRSlope = BR.Slope
) %>%
mutate(
dateTime = as.POSIXct(ymd_hms(dateTime)),
Time = as_hms(ymd_hms(dateTime)),
Date = as.Date(dateTime, format = "%Y/%m/%d"),
)
tb_respirometry <- tb_respirometry %>%
mutate(
volume_net = volume_ch - mass,
MR_wBR = abs(slope_wBR)*(volume_net/1000)*60*60, #uncomment for mgO2/hr instead
BR = BRSlope*(volume_ch/1000)*60*60,
MR = MR_wBR + BR,
BR_perc = (BR/MR_wBR)*100
)
#Calculate RMR and variance in MR for each fish, create a new table for individual RMR calcs
tb_rmr <- tb_respirometry %>%
group_by(chamber_ch) %>%
arrange(chamber_ch, dateTime) %>%
slice(3:(n() - 1)) %>%
ungroup() %>%
group_by(
ID_fish, mass, length, volume_net, chamber_ch, Temp_class
) %>%
arrange(ID_fish) %>%
drop_na() %>%
summarise(
RMR = mean(MR),
RMR_var = var(MR),
RMR_perc = (RMR_var/RMR)*100,
time_start = dateTime %>% min(),
time_end = dateTime %>% max()
)
#Calculate SMR for each fish, create a new table for individual SMR calcs
tb_smr <-
tb_respirometry %>%
group_by(
ID_fish, mass,length, volume_net, chamber_ch
) %>%
arrange(ID_fish) %>%
drop_na() %>%
summarise(
SMR = calcSMR(MR)$low10pc %>% unname(),
time_start = dateTime %>% min(),
time_end = dateTime %>% max()
)
# Process XLSX data
tb_mmr <- read_excel(current_files$xlsx) %>%
rename(
chamber_ch = Chamber.No,
ID_fish = Ind,
mass = Mass,
volume_ch = Ch.Volume,
DOunit = DO.unit,
dateTime = Date.Time,
Temp_class = Temp.class,
slope_wBR = Slope.with.BR,
BRSlope = BR.Slope
) %>%
mutate(
dateTime = as.POSIXct(ymd_hms(dateTime)),
Time = as_hms(ymd_hms(dateTime)),
Date = as.Date(dateTime, format = "%Y/%m/%d"),
) %>%
drop_na()
tb_mmr <- tb_mmr %>%
mutate(
volume_net = volume_ch - mass,
MMR_wBR = abs(slope_wBR)*(volume_net/1000)*60*60, #uncomment for mgO2/hr instead
BR = abs(BRSlope)*(volume_ch/1000)*60*60,
MMR = MMR_wBR - BR,
BR_perc = (BR/MMR_wBR)*100
) %>%
arrange(ID_fish)
# Combine the processed data
# Combine the processed data, selecting only specific columns
tb_MR_master <- tb_rmr %>%
left_join(select(tb_smr, ID_fish, SMR), by = "ID_fish") %>%
left_join(select(tb_mmr, ID_fish, MMR), by = "ID_fish") %>%
select(-matches(".x$")) %>%
rename_with(~gsub("\\.y$", "", .), matches(".y$"))
# Export data
setwd(MRcalcs_wd)
write.csv(tb_MR_master, file = paste0(date, "_MRcalcs.csv"), col.names = NA, row.names = FALSE)
}
# read all "_MRcalcs.csv" files, bind into one table, and add columns with logged values
file.list <- list.files(MRcalcs_wd)
setwd(MRcalcs_wd)
calcs_all <- lapply(file.list, read_csv)
tb_MRcalcs <- bind_rows(calcs_all, .id="Resp_Day") %>%
mutate(
log_SMR_low10pc = log10(SMR),
log_mass = log10(mass),
log_RMR = log10(RMR),
log_MMR = log10(MMR),
Resp_Day = as.numeric(Resp_Day)
)
# assign Month category to all timepoints
tb_MRcalcs$Month <- cut(
tb_MRcalcs$Resp_Day,
breaks = c(0, 13, seq(23, max(tb_MRcalcs$Resp_Day) + 10, by = 10)),
right = FALSE,
labels = c("April", "May", "June", "July", "August", "September")
)
setwd(raw_data_wd)
biometrics <- read.csv("FinalBiometrics.csv")
ds <- left_join(tb_MRcalcs, biometrics, by = c("ID_fish")) #a table with all resp data and metadata collected in dissection
ds1 <- ds[which(ds$Month != "April"),] # remove April data from analysis due to temperature inconsistencyMultivariate mixed-effects models to test the relationships between growth and metabolism in Galaxias maculatus.
ds1 <- ds1 %>%
# Filter to keep only ID_fish where a mass exists for Month == "May"
filter(ID_fish %in% ID_fish[Month == "May" & !is.na(mass)])
month_order <- c("May", "June", "July", "August", "September") # assign chronological order to Month (otherwise alphabetical)
ds1$Month <- factor(ds1$Month, levels = month_order) # assign month_order to factor variable "Month" in ds1
month_order_pop <- c("April", "May", "June", "July", "August", "September") # assign chronological order to Month (otherwise alphabetical)
min_mass = min(ds1$mass)
min_log_mass = min(ds1$log_mass)
min_finalmass = min(((ds1[which(!is.na(ds1$Final_Mass_g)),]$Final_Mass_g)))
min_log_finalmass = min((log10(ds1[which(!is.na(ds1$Final_Mass_g)),]$Final_Mass_g)))
# format data columns and generate centred and log-scaled variables
ds1 <- ds1 %>%
mutate(
mass_centre = mass - min_mass, # centred mass
log_mass_centre = log_mass - log10(min_mass), # log-scaled centred mass
Temp_class = as.factor(Temp_class), # set temperature treatments as factors
SMR1 = log10(SMR), # log-scaled SMR
RMR1 = log10(RMR), # log-scaled RMR
MMR1 = log10(MMR), # log-scaled MMR
mass1 = log_mass, # log-scaled mass
Month_int = case_when(Month == "May" ~ 1,
Month == "June" ~ 2,
Month == "July" ~ 3,
Month == "August" ~ 4,
Month == "September" ~ 5),
Month_centred = Month_int - 1, # left-centred mass so May is the intercept
Month_centred_factor = as.factor(Month_centred), # set centred month as factors
Rack = as.integer(Rack), # set rack ID to integer (for now)
Bin = as.integer(Bin), # set bin ID to integer (for now)
finalmass_centre = Final_Mass_g - min_finalmass, # centre last mass measurement taken
log_finalmass_centre = log10(Final_Mass_g) - min_log_finalmass, # log-scaled and centred final mass
log_fatmass = log10(Fat_mass_g) # log-scaled fat mass
)
ds1 <- tidyr::unite(data = ds1, Rack, Bin, col = "tankID", sep = "_") # merge rack and bin data into one identifier, call it tankID
#read in density data table, call it BinCounts
setwd(raw_data_wd)
BinCounts <- read_excel("BinCounts.xlsx", sheet = "data")
tb_density <- BinCounts %>%
group_by_(
"Month",
"Rack") %>%
summarise(
Population = sum(Population_tagged)) %>%
mutate(
Month = factor(Month, levels = month_order_pop))
ggplot(BinCounts, aes(x = factor(Month, levels = month_order_pop), y = Population_tagged, group = tankID, color = factor(tankID))) +
geom_line() +
geom_point() +
geom_jitter(height = 0.1, width = 0.25) +
facet_grid(.~Rack) +
theme_classic()
#merge BinCounts and ds1
ds1 <- left_join(ds1, BinCounts %>% dplyr::select("Month", "tankID", "Population"), by = c("Month", "tankID"))
#reformat Sex variable so that both unknown and immature fish are "NA"
ds1$Sex <- ifelse(ds1$Sex == "I", NA, ds1$Sex)
# rename "population" as density, set tankID as factor
ds1 <- ds1 %>%
rename(
Density = Population
) %>%
mutate(
tankID = as.factor(tankID)
)
# check structure of ds1
str(ds1)## tibble [624 × 41] (S3: tbl_df/tbl/data.frame)
## $ Resp_Day : num [1:624] 13 13 13 13 13 13 13 13 13 13 ...
## $ ID_fish : num [1:624] 81 82 83 84 85 86 87 88 89 90 ...
## $ chamber_ch : chr [1:624] "A1" "B3" "B4" "A2" ...
## $ Temp_class : Factor w/ 2 levels "18","23": 2 2 2 2 2 2 2 2 2 2 ...
## $ RMR : num [1:624] 0.952 0.412 0.712 1.016 0.876 ...
## $ RMR_var : num [1:624] 0.0149 0.0135 0.0429 0.0696 0.0454 ...
## $ RMR_perc : num [1:624] 1.57 3.28 6.02 6.85 5.18 ...
## $ time_start : POSIXct[1:624], format: "2023-05-22 17:53:00" "2023-05-22 17:53:00" ...
## $ time_end : POSIXct[1:624], format: "2023-05-23 08:03:00" "2023-05-23 08:03:00" ...
## $ mass : num [1:624] 2.65 0.81 1.78 2.73 2.29 0.98 2.04 2.43 2.31 1.65 ...
## $ length : num [1:624] 77 62 73 78 80 64 80 77 78 68 ...
## $ volume_net : num [1:624] 707 709 708 707 708 ...
## $ SMR : num [1:624] 0.845 0.284 0.53 0.804 0.707 ...
## $ MMR : num [1:624] 1.43 1.11 1.44 1.64 1.69 ...
## $ log_SMR_low10pc : num [1:624] -0.0733 -0.5461 -0.2758 -0.0948 -0.1506 ...
## $ log_mass : num [1:624] 0.4232 -0.0915 0.2504 0.4362 0.3598 ...
## $ log_RMR : num [1:624] -0.02124 -0.38552 -0.14722 0.00674 -0.05726 ...
## $ log_MMR : num [1:624] 0.1544 0.0458 0.1591 0.2154 0.2275 ...
## $ Month : chr [1:624] "May" "May" "May" "May" ...
## $ tankID : Factor w/ 20 levels "1_10","1_11",..: 17 17 17 17 17 17 17 17 18 18 ...
## $ TagID : chr [1:624] "YL" "YR" "BrL" "BrR" ...
## $ Final_Mass_g : num [1:624] 4.62 3.65 2.64 6.33 5.19 ...
## $ Final_Length_mm : int [1:624] 94 86 84 98 100 95 90 89 91 98 ...
## $ Gonad_mass_g : num [1:624] 0.035 0.022 0.018 0.039 0.052 0.028 0.362 0.033 0.015 0.093 ...
## $ Sex : chr [1:624] "F" "M" "F" "M" ...
## $ Fat_mass_g : num [1:624] 0.315 0.123 0.062 0.314 0.214 0.091 0.076 0.044 0.15 0.209 ...
## $ Ventricle_mass_g : chr [1:624] "0.017" "" "" "" ...
## $ Otolith_age_y : num [1:624] 1 NA NA NA NA 1 1 NA NA 1 ...
## $ mass_centre : num [1:624] 2.07 0.23 1.2 2.15 1.71 0.4 1.46 1.85 1.73 1.07 ...
## $ log_mass_centre : num [1:624] 0.66 0.145 0.487 0.673 0.596 ...
## $ SMR1 : num [1:624] -0.0733 -0.5461 -0.2758 -0.0948 -0.1506 ...
## $ RMR1 : num [1:624] -0.02124 -0.38552 -0.14722 0.00674 -0.05726 ...
## $ MMR1 : num [1:624] 0.1544 0.0458 0.1591 0.2154 0.2275 ...
## $ mass1 : num [1:624] 0.4232 -0.0915 0.2504 0.4362 0.3598 ...
## $ Month_int : num [1:624] 1 1 1 1 1 1 1 1 1 1 ...
## $ Month_centred : num [1:624] 0 0 0 0 0 0 0 0 0 0 ...
## $ Month_centred_factor: Factor w/ 5 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ finalmass_centre : num [1:624] 3.3 2.33 1.32 5.01 3.87 ...
## $ log_finalmass_centre: num [1:624] 0.544 0.442 0.301 0.681 0.594 ...
## $ log_fatmass : num [1:624] -0.502 -0.91 -1.208 -0.503 -0.67 ...
## $ Density : num [1:624] 8 8 8 8 8 8 8 8 5 5 ...head(ds1)## # A tibble: 6 × 41
## Resp_Day ID_fish chamber_ch Temp_class RMR RMR_var RMR_perc
## <dbl> <dbl> <chr> <fct> <dbl> <dbl> <dbl>
## 1 13 81 A1 23 0.952 0.0149 1.57
## 2 13 82 B3 23 0.412 0.0135 3.28
## 3 13 83 B4 23 0.712 0.0429 6.02
## 4 13 84 A2 23 1.02 0.0696 6.85
## 5 13 85 A4 23 0.876 0.0454 5.18
## 6 13 86 A3 23 0.751 0.0504 6.71
## # ℹ 34 more variables: time_start <dttm>, time_end <dttm>, mass <dbl>,
## # length <dbl>, volume_net <dbl>, SMR <dbl>, MMR <dbl>,
## # log_SMR_low10pc <dbl>, log_mass <dbl>, log_RMR <dbl>, log_MMR <dbl>,
## # Month <chr>, tankID <fct>, TagID <chr>, Final_Mass_g <dbl>,
## # Final_Length_mm <int>, Gonad_mass_g <dbl>, Sex <chr>, Fat_mass_g <dbl>,
## # Ventricle_mass_g <chr>, Otolith_age_y <dbl>, mass_centre <dbl>,
## # log_mass_centre <dbl>, SMR1 <dbl>, RMR1 <dbl>, MMR1 <dbl>, mass1 <dbl>, …# check sex ratio of population
tb_sexratio <- ds1[which(ds1$Month == "September"),] %>%
group_by_(
"Temp_class",
"Sex"
) %>%
count()
sexorder <- c("Immature", "Female", "Male")
tb_sexratio$Sex <- as.character(tb_sexratio$Sex)
tb_sexratio$Sex[tb_sexratio$Sex==""] <- NA
tb_sexratio$Sex[is.na(tb_sexratio$Sex)] <- "Immature"
tb_sexratio$Sex[tb_sexratio$Sex == "F"] <- "Female"
tb_sexratio$Sex[tb_sexratio$Sex == "M"] <- "Male"
tb_sexratio <- tb_sexratio %>%
mutate(
Sex = factor(Sex, levels = sexorder))
fig3b1 <- ggplot(tb_sexratio, aes(x = " ", y = n, fill = Sex)) +
geom_bar(position = "stack", stat = "identity") +
scale_fill_manual(values = c("gold3", "darkgreen", "purple")) +
labs(x = "Whole Population", y = "Count (absolute)", fill = "Sex") +
theme_classic()
fig3b2 <- ggplot(tb_sexratio, aes(x = factor(Temp_class), y = n, fill = Sex)) +
geom_bar(position = "fill", stat = "identity") +
scale_fill_manual(values = c("gold3", "darkgreen", "purple")) +
labs(x = expression("Temperature ("*degree*C*")"), y = "Proportion", fill = "Sex") +
theme_classic()
fig3b <- ggarrange(print(fig3b1), print(fig3b2), ncol = 2, nrow = 1, widths = c(1, 2), align = "v", common.legend = TRUE, legend = "right")
fig3b
GROWTH ANALYSES
Visualize growth data
hist(ds1$mass) # plot all mass values to see mass range and most frequent values
ds1$Month <- factor(ds1$Month, levels = month_order) # reassign month_order to factor variable "Month" in ds1
ds1$Sex[ds1$Sex == ""] <- NA
ds1$Sex <- factor(ds1$Sex)
# Plot mass over time, connect by ID_fish to show individual trends
ggplot(ds1, aes(x = Month, y = mass, group = ID_fish)) + # plot raw data by ID
geom_point(size = 2) + # size of data points
geom_line() + # connect points by individual ID
facet_wrap(. ~ Temp_class) + # produces one plot per temperature treatment
theme_bw()
Log-linearize mass data
We can also visualise this by plotting log10 mass over log10_x (note that we cannot use log10_y as some masses were < 1 initially)
ggplot(ds1, aes(x = Month_int, y = log_mass, group = ID_fish)) + # plot logged mass data by ID
geom_point(size = 2) + # size of data points
geom_line() + # connect points by individual ID
facet_wrap(. ~ Temp_class) + # produces one plot per temperature treatment
theme_bw() +
scale_x_log10() # display x-axis on log scale
Univariate growth models for two temperature treatments
ds18 <- ds1[which(ds1$Temp_class == "18"),] # subset ds for only fish at 18C
ds23 <- ds1[which(ds1$Temp_class == "23"),] # subset ds for only fish at 23C
mass18_mod<-lme4::lmer(log_mass ~ Month + (1 + Month_int|ID_fish), # model structure, 18C only
data=ds18, REML = T)
summary(mass18_mod) # generate model output## Linear mixed model fit by REML ['lmerMod']
## Formula: log_mass ~ Month + (1 + Month_int | ID_fish)
## Data: ds18
##
## REML criterion at convergence: -492
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8214 -0.4648 0.0139 0.4865 2.3450
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.040770 0.20191
## Month_int 0.002038 0.04515 -0.57
## Residual 0.003528 0.05939
## Number of obs: 317, groups: ID_fish, 72
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.32789 0.02232 14.689
## MonthJune 0.15461 0.01147 13.478
## MonthJuly 0.22113 0.01498 14.761
## MonthAugust 0.15515 0.01989 7.799
## MonthSeptember 0.36370 0.02487 14.622
##
## Correlation of Fixed Effects:
## (Intr) MnthJn MnthJl MnthAg
## MonthJune -0.364
## MonthJuly -0.412 0.629
## MonthAugust -0.410 0.600 0.770
## MonthSptmbr -0.408 0.581 0.776 0.865plot(mass18_mod) # residuals plot
mass23_mod<-lme4::lmer(log_mass ~ Month + (1 + Month_int|ID_fish), # model structure, 23C only
data=ds23, REML=T)
summary(mass23_mod) # generate model output## Linear mixed model fit by REML ['lmerMod']
## Formula: log_mass ~ Month + (1 + Month_int | ID_fish)
## Data: ds23
##
## REML criterion at convergence: -515.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9289 -0.3941 -0.0018 0.5057 4.8187
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.040995 0.20247
## Month_int 0.001445 0.03802 -0.46
## Residual 0.003138 0.05602
## Number of obs: 307, groups: ID_fish, 64
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.21996 0.02451 8.974
## MonthJune 0.12091 0.01123 10.771
## MonthJuly 0.21448 0.01402 15.297
## MonthAugust 0.22778 0.01773 12.848
## MonthSeptember 0.35157 0.02193 16.028
##
## Correlation of Fixed Effects:
## (Intr) MnthJn MnthJl MnthAg
## MonthJune -0.298
## MonthJuly -0.335 0.623
## MonthAugust -0.341 0.611 0.773
## MonthSptmbr -0.337 0.589 0.778 0.858plot(mass23_mod) # residuals plot
Combined growth model
# mass model including both temperature treatments, and their interaction with Month
mass_mod <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),
data=ds1, REML=T)
summary(mass_mod) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (1 +
## Month_int | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -732.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6965 -0.4297 0.0123 0.4659 4.5929
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.032801 0.18111
## Month_int 0.001623 0.04029 -0.54
## Residual 0.003489 0.05906
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.287001 0.051703 379.511550 5.551
## MonthJune 0.157704 0.012697 354.616467 12.420
## MonthJuly 0.228444 0.016168 237.622029 14.130
## MonthAugust 0.167240 0.021107 155.507265 7.923
## MonthSeptember 0.379389 0.025536 111.141150 14.857
## as.factor(Temp_class)23 -0.059851 0.037534 99.223469 -1.595
## SexM -0.073375 0.033162 88.017749 -2.213
## Density 0.011296 0.005978 429.998459 1.890
## MonthJune:as.factor(Temp_class)23 -0.024598 0.019947 352.647746 -1.233
## MonthJuly:as.factor(Temp_class)23 -0.010083 0.025097 226.080174 -0.402
## MonthAugust:as.factor(Temp_class)23 0.060749 0.032152 137.635190 1.889
## MonthSeptember:as.factor(Temp_class)23 -0.023724 0.039367 100.653936 -0.603
## Pr(>|t|)
## (Intercept) 5.35e-08 ***
## MonthJune < 2e-16 ***
## MonthJuly < 2e-16 ***
## MonthAugust 4.15e-13 ***
## MonthSeptember < 2e-16 ***
## as.factor(Temp_class)23 0.1140
## SexM 0.0295 *
## Density 0.0595 .
## MonthJune:as.factor(Temp_class)23 0.2183
## MonthJuly:as.factor(Temp_class)23 0.6882
## MonthAugust:as.factor(Temp_class)23 0.0609 .
## MonthSeptember:as.factor(Temp_class)23 0.5481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MnthJn MnthJl MnthAg MnthSp a.(T_) SexM Densty
## MonthJune -0.281
## MonthJuly -0.384 0.632
## MonthAugust -0.473 0.615 0.787
## MonthSptmbr -0.423 0.599 0.794 0.873
## as.fc(T_)23 -0.358 0.239 0.273 0.282 0.278
## SexM -0.215 -0.004 -0.006 -0.007 -0.006 -0.098
## Density -0.861 0.139 0.239 0.345 0.285 0.116 -0.021
## MnthJn:.(T_)23 0.143 -0.631 -0.392 -0.377 -0.370 -0.364 0.002 -0.047
## MnthJl:.(T_)23 0.219 -0.402 -0.636 -0.496 -0.502 -0.414 0.003 -0.121
## MnA:.(T_)23 0.288 -0.400 -0.511 -0.648 -0.566 -0.425 0.004 -0.201
## MnS:.(T_)23 0.245 -0.384 -0.507 -0.555 -0.639 -0.419 0.003 -0.151
## MnthJn:.(T_)23 MnthJl:.(T_)23 MA:.(T
## MonthJune
## MonthJuly
## MonthAugust
## MonthSptmbr
## as.fc(T_)23
## SexM
## Density
## MnthJn:.(T_)23
## MnthJl:.(T_)23 0.625
## MnA:.(T_)23 0.608 0.778
## MnS:.(T_)23 0.591 0.785 0.864plot(mass_mod) # residuals plot
ranova(mass_mod) # ANOVA table output to test random effects## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + Month:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 16 366.13 -700.25
## Month_int in (1 + Month_int | ID_fish) 14 305.38 -582.77 121.48 2 < 2.2e-16
##
## <none>
## Month_int in (1 + Month_int | ID_fish) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# mass model including both temperature treatments, and their interaction with Month
mass_mod_noMonth <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (1|ID_fish),
data=ds1, REML=T)
summary(mass_mod_noMonth) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (1 |
## ID_fish)
## Data: ds1
##
## REML criterion at convergence: -610.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6634 -0.5041 0.0313 0.4740 3.9635
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.022746 0.15082
## Residual 0.007567 0.08699
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.304207 0.053767 382.535854 5.658
## MonthJune 0.157025 0.016791 355.224021 9.352
## MonthJuly 0.227033 0.017176 361.268479 13.218
## MonthAugust 0.164380 0.018271 377.508492 8.997
## MonthSeptember 0.376529 0.018271 377.508492 20.608
## as.factor(Temp_class)23 -0.061639 0.037779 137.305730 -1.632
## SexM -0.072332 0.033449 88.014836 -2.162
## Density 0.008948 0.006287 440.998220 1.423
## MonthJune:as.factor(Temp_class)23 -0.024244 0.026455 353.637257 -0.916
## MonthJuly:as.factor(Temp_class)23 -0.008998 0.026629 355.432235 -0.338
## MonthAugust:as.factor(Temp_class)23 0.063283 0.027233 362.066784 2.324
## MonthSeptember:as.factor(Temp_class)23 -0.021386 0.027107 360.854223 -0.789
## Pr(>|t|)
## (Intercept) 3.01e-08 ***
## MonthJune < 2e-16 ***
## MonthJuly < 2e-16 ***
## MonthAugust < 2e-16 ***
## MonthSeptember < 2e-16 ***
## as.factor(Temp_class)23 0.1051
## SexM 0.0333 *
## Density 0.1554
## MonthJune:as.factor(Temp_class)23 0.3601
## MonthJuly:as.factor(Temp_class)23 0.7356
## MonthAugust:as.factor(Temp_class)23 0.0207 *
## MonthSeptember:as.factor(Temp_class)23 0.4307
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MnthJn MnthJl MnthAg MnthSp a.(T_) SexM Densty
## MonthJune -0.248
## MonthJuly -0.355 0.506
## MonthAugust -0.505 0.495 0.538
## MonthSptmbr -0.505 0.495 0.538 0.588
## as.fc(T_)23 -0.355 0.230 0.241 0.250 0.250
## SexM -0.206 -0.004 -0.007 -0.010 -0.010 -0.099
## Density -0.871 0.110 0.237 0.419 0.419 0.121 -0.023
## MnthJn:.(T_)23 0.129 -0.631 -0.313 -0.300 -0.300 -0.353 0.002 -0.037
## MnthJl:.(T_)23 0.200 -0.323 -0.637 -0.333 -0.333 -0.360 0.004 -0.120
## MnA:.(T_)23 0.311 -0.328 -0.353 -0.657 -0.381 -0.368 0.006 -0.249
## MnS:.(T_)23 0.296 -0.328 -0.350 -0.375 -0.652 -0.368 0.005 -0.231
## MnthJn:.(T_)23 MnthJl:.(T_)23 MA:.(T
## MonthJune
## MonthJuly
## MonthAugust
## MonthSptmbr
## as.fc(T_)23
## SexM
## Density
## MnthJn:.(T_)23
## MnthJl:.(T_)23 0.499
## MnA:.(T_)23 0.492 0.509
## MnS:.(T_)23 0.494 0.510 0.529plot(mass_mod_noMonth) # residuals plot
ranova(mass_mod_noMonth) # ANOVA table output to test random effects## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + Month:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 14 305.38 -582.77
## (1 | ID_fish) 13 121.89 -217.78 366.99 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# mass model with random intercepts by temp treatment
mass_mod_trtint <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (0 + Month_int|Temp_class),
data=ds1, REML=T)
summary(mass_mod_trtint) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (0 +
## Month_int | Temp_class)
## Data: ds1
##
## REML criterion at convergence: -243.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.86639 -0.70696 -0.01277 0.67903 2.94871
##
## Random effects:
## Groups Name Variance Std.Dev.
## Temp_class Month_int 0.03024 0.1739
## Residual 0.03013 0.1736
## Number of obs: 453, groups: Temp_class, 2
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.301995 0.181676 441.000000 1.662
## MonthJune 0.159454 0.177058 441.000000 0.901
## MonthJuly 0.229581 0.349405 441.000000 0.657
## MonthAugust 0.164793 0.522797 441.000000 0.315
## MonthSeptember 0.376942 0.696419 441.000000 0.541
## as.factor(Temp_class)23 -0.061299 0.248776 441.000000 -0.246
## SexM -0.073232 0.016703 441.000000 -4.384
## Density 0.009287 0.006273 441.000000 1.480
## MonthJune:as.factor(Temp_class)23 -0.026626 0.251518 441.000000 -0.106
## MonthJuly:as.factor(Temp_class)23 -0.011500 0.494684 441.000000 -0.023
## MonthAugust:as.factor(Temp_class)23 0.062917 0.739690 441.000000 0.085
## MonthSeptember:as.factor(Temp_class)23 -0.021724 0.985139 441.000000 -0.022
## Pr(>|t|)
## (Intercept) 0.0972 .
## MonthJune 0.3683
## MonthJuly 0.5115
## MonthAugust 0.7527
## MonthSeptember 0.5886
## as.factor(Temp_class)23 0.8055
## SexM 1.46e-05 ***
## Density 0.1395
## MonthJune:as.factor(Temp_class)23 0.9157
## MonthJuly:as.factor(Temp_class)23 0.9815
## MonthAugust:as.factor(Temp_class)23 0.9323
## MonthSeptember:as.factor(Temp_class)23 0.9824
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MnthJn MnthJl MnthAg MnthSp a.(T_) SexM Densty
## MonthJune 0.920
## MonthJuly 0.941 0.987
## MonthAugust 0.946 0.986 0.996
## MonthSptmbr 0.949 0.986 0.997 0.998
## as.fc(T_)23 -0.686 -0.674 -0.689 -0.693 -0.695
## SexM -0.021 -0.001 -0.001 -0.001 -0.001 -0.008
## Density -0.257 0.011 0.012 0.015 0.011 0.018 -0.047
## MnthJn:.(T_)23 -0.649 -0.704 -0.694 -0.694 -0.694 0.944 0.001 -0.004
## MnthJl:.(T_)23 -0.665 -0.697 -0.706 -0.704 -0.704 0.972 0.001 -0.007
## MnA:.(T_)23 -0.669 -0.697 -0.704 -0.707 -0.706 0.978 0.000 -0.009
## MnS:.(T_)23 -0.671 -0.697 -0.704 -0.706 -0.707 0.981 0.000 -0.006
## MnthJn:.(T_)23 MnthJl:.(T_)23 MA:.(T
## MonthJune
## MonthJuly
## MonthAugust
## MonthSptmbr
## as.fc(T_)23
## SexM
## Density
## MnthJn:.(T_)23
## MnthJl:.(T_)23 0.983
## MnA:.(T_)23 0.983 0.996
## MnS:.(T_)23 0.982 0.996 0.998plot(mass_mod_trtint) # residuals plot
ranova(mass_mod_trtint) # ANOVA table output to test random effects## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (0 + Month_int | Temp_class) + Month:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 14 121.89 -215.78
## Month_int in (0 + Month_int | Temp_class) 14 121.89 -215.78 -1.7053e-13 0
## Pr(>Chisq)
## <none>
## Month_int in (0 + Month_int | Temp_class)mass_mod_sex <- lmer(log_mass ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):Month + Density + (1 + Month_int|ID_fish),
data=ds1, REML=T)
summary(mass_mod_sex)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log_mass ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):Month +
## Density + (1 + Month_int | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -706.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4000 -0.4181 0.0150 0.4378 4.4184
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.033320 0.18254
## Month_int 0.001501 0.03874 -0.56
## Residual 0.003432 0.05858
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error
## Density 1.228e-02 5.915e-03
## as.factor(Sex)F:as.factor(Temp_class)18 2.817e-01 5.316e-02
## as.factor(Sex)M:as.factor(Temp_class)18 2.026e-01 5.901e-02
## as.factor(Sex)F:as.factor(Temp_class)23 2.237e-01 5.597e-02
## as.factor(Sex)M:as.factor(Temp_class)23 1.435e-01 5.888e-02
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJune 1.598e-01 1.571e-02
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJune 1.544e-01 2.050e-02
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJune 1.078e-01 2.098e-02
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJune 1.617e-01 2.225e-02
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJuly 2.577e-01 1.971e-02
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJuly 1.797e-01 2.557e-02
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJuly 1.689e-01 2.602e-02
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJuly 2.739e-01 2.757e-02
## as.factor(Sex)F:as.factor(Temp_class)18:MonthAugust 1.978e-01 2.530e-02
## as.factor(Sex)M:as.factor(Temp_class)18:MonthAugust 1.170e-01 3.257e-02
## as.factor(Sex)F:as.factor(Temp_class)23:MonthAugust 1.757e-01 3.274e-02
## as.factor(Sex)M:as.factor(Temp_class)23:MonthAugust 2.868e-01 3.466e-02
## as.factor(Sex)F:as.factor(Temp_class)18:MonthSeptember 4.041e-01 3.066e-02
## as.factor(Sex)M:as.factor(Temp_class)18:MonthSeptember 3.394e-01 3.983e-02
## as.factor(Sex)F:as.factor(Temp_class)23:MonthSeptember 2.912e-01 4.032e-02
## as.factor(Sex)M:as.factor(Temp_class)23:MonthSeptember 4.282e-01 4.265e-02
## df t value
## Density 4.209e+02 2.076
## as.factor(Sex)F:as.factor(Temp_class)18 3.488e+02 5.299
## as.factor(Sex)M:as.factor(Temp_class)18 2.778e+02 3.434
## as.factor(Sex)F:as.factor(Temp_class)23 2.403e+02 3.997
## as.factor(Sex)M:as.factor(Temp_class)23 2.379e+02 2.437
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJune 3.451e+02 10.173
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJune 3.442e+02 7.530
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJune 3.436e+02 5.136
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJune 3.443e+02 7.268
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJuly 2.361e+02 13.070
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJuly 2.266e+02 7.028
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJuly 2.229e+02 6.491
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJuly 2.238e+02 9.937
## as.factor(Sex)F:as.factor(Temp_class)18:MonthAugust 1.475e+02 7.821
## as.factor(Sex)M:as.factor(Temp_class)18:MonthAugust 1.366e+02 3.593
## as.factor(Sex)F:as.factor(Temp_class)23:MonthAugust 1.287e+02 5.365
## as.factor(Sex)M:as.factor(Temp_class)23:MonthAugust 1.292e+02 8.273
## as.factor(Sex)F:as.factor(Temp_class)18:MonthSeptember 1.054e+02 13.181
## as.factor(Sex)M:as.factor(Temp_class)18:MonthSeptember 9.942e+01 8.523
## as.factor(Sex)F:as.factor(Temp_class)23:MonthSeptember 9.524e+01 7.221
## as.factor(Sex)M:as.factor(Temp_class)23:MonthSeptember 9.541e+01 10.040
## Pr(>|t|)
## Density 0.038492 *
## as.factor(Sex)F:as.factor(Temp_class)18 2.07e-07 ***
## as.factor(Sex)M:as.factor(Temp_class)18 0.000686 ***
## as.factor(Sex)F:as.factor(Temp_class)23 8.53e-05 ***
## as.factor(Sex)M:as.factor(Temp_class)23 0.015541 *
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJune < 2e-16 ***
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJune 4.50e-13 ***
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJune 4.70e-07 ***
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJune 2.46e-12 ***
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJuly < 2e-16 ***
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJuly 2.43e-11 ***
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJuly 5.45e-10 ***
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJuly < 2e-16 ***
## as.factor(Sex)F:as.factor(Temp_class)18:MonthAugust 9.26e-13 ***
## as.factor(Sex)M:as.factor(Temp_class)18:MonthAugust 0.000455 ***
## as.factor(Sex)F:as.factor(Temp_class)23:MonthAugust 3.63e-07 ***
## as.factor(Sex)M:as.factor(Temp_class)23:MonthAugust 1.39e-13 ***
## as.factor(Sex)F:as.factor(Temp_class)18:MonthSeptember < 2e-16 ***
## as.factor(Sex)M:as.factor(Temp_class)18:MonthSeptember 1.75e-13 ***
## as.factor(Sex)F:as.factor(Temp_class)23:MonthSeptember 1.26e-10 ***
## as.factor(Sex)M:as.factor(Temp_class)23:MonthSeptember < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(mass_mod_sex)
Plot predicted mean level trend in mass over time, by temperature
pred <- ggpredict(mass_mod, terms = c("Month", "Temp_class")) # generate trend lines for the two temperature treatments
ggplot() +
geom_errorbar(
data = pred,
aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) + # add 95% confidence interval bars
geom_point(
data = pred,
aes(x = x, y = predicted, color = group, group = group)) + # predicted level by temp*month
geom_line(
data = pred,
aes(x = x, y = predicted, color = group, group = group)) + # show level trend over time
labs(
x = "Month",
y = "Predicted log_mass",
color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00"))
F test on interaction effect
anova(mass_mod, type = "3") # type 3 anova because we expect a meaningful interaction## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Month 1.68541 0.42135 4 205.57 120.7818 < 2.2e-16 ***
## as.factor(Temp_class) 0.01090 0.01090 1 88.16 3.1235 0.08063 .
## Sex 0.01708 0.01708 1 88.02 4.8958 0.02951 *
## Density 0.01246 0.01246 1 430.00 3.5706 0.05948 .
## Month:as.factor(Temp_class) 0.10282 0.02571 4 201.73 7.3687 1.472e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Generate fitted data
z_mass<-augment(mass_mod) # get variables and fitted values from model
head(z_mass)## # A tibble: 6 × 19
## .rownames log_mass Month `as.factor(Temp_class)` Sex Density Month_int
## <chr> <dbl> <fct> <fct> <fct> <dbl> <dbl>
## 1 1 0.423 May 23 F 8 1
## 2 2 -0.0915 May 23 M 8 1
## 3 3 0.250 May 23 F 8 1
## 4 4 0.436 May 23 M 8 1
## 5 5 0.360 May 23 M 8 1
## 6 7 0.310 May 23 M 8 1
## # ℹ 12 more variables: ID_fish <dbl>, .fitted <dbl>, .resid <dbl>, .hat <dbl>,
## # .cooksd <dbl>, .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>,
## # .sqrtrwt <dbl>, .weights <dbl>, .wtres <dbl>##### check assumptions: plot all residual values against the fitted
ggplot(z_mass, aes(x = .fitted, y = .resid)) + # same plot as before
geom_point(size = 2) 
FIGURE 1
fig1.1 <- ggplot(z_mass, aes(x = Month, y = log_mass, group = ID_fish)) +
geom_line(aes(y = .fitted, color = `as.factor(Temp_class)`), # note this will look odd if there are additional fixed effects
alpha = 0.6) +
theme_classic() +
labs(x = "Month", y = expression(log[10]~(mass~(g))), color = expression("Temperature ("*degree*C*")")) +
ylim(-0.24,1) +
scale_color_manual(values = c("#78B7C5", "#F21A00"))
fig1.2 <- ggplot() +
geom_errorbar(
data = pred,
aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2, alpha = 0.3) + # add 95% confidence interval bars
geom_line(
data = pred,
aes(x = x, y = predicted, color = group, group = group),
size = 1) +
labs(
x = "Month",
y = expression(log[10]~(mass~(g))),
color = expression("Temperature ("*degree*C*")")) +
theme_classic() +
ylim(-0.24,1) +
scale_color_manual(values = c("#78B7C5", "#F21A00"))
fig1 <- ggarrange(print(fig1.1 + rremove("xlab") + rremove ("ylab")),
print(fig1.2 + rremove("xlab") + rremove ("ylab")),
ncol = 2, nrow = 1, widths = c(1, 1.1),
align = "v",
common.legend = TRUE,
legend = "right")
require(grid)
jpeg(file = "Fig1.jpeg", width=960, height=480, units = "px")
annotate_figure(
fig1,
bottom = textGrob("Month"),
left = textGrob(
expression(log[10]~(mass~(g))),
rot = 90, # Rotate the y-axis label 90 degrees
gp = gpar(fontsize = 12) # Optionally adjust the font size
)
)
setwd(figures_wd)
jpeg(file = "Fig1.jpeg", width=960, height=480, units = "px")
annotate_figure(
fig1,
bottom = textGrob("Month"),
left = textGrob(
expression(log[10]~(mass~(g))),
rot = 90, # Rotate the y-axis label 90 degrees
gp = gpar(fontsize = 12) # Optionally adjust the font size
)
)
dev.off()## png
## 2SMR MODEL
Visualize SMR data
hist(ds1$SMR) # plot all raw SMR values in a histogram, note skewness
hist(log10(ds1$SMR)) # plot all log-SMR values in a histogram
ggplot(ds1, aes(x = mass, y = SMR, group = ID_fish)) + # plot raw SMR data by ID against mass
geom_point(size = 2) + # size of data points
facet_wrap(. ~ Temp_class) + # produces one plot per temperature treatment
theme_bw() 
Visualize log-mass and SMR relationship
ggplot(ds1, aes(x = log_mass, y = SMR, group = ID_fish)) + # plot raw SMR data by ID against log10(mass)
geom_point(size = 2) + # size of data points
facet_wrap(. ~ Temp_class) + # produces one plot per temperature treatment
theme_bw() 
Note that this relationship is non-linear, we expect this to be described by a power function y = a*x^b
This reflects our understanding of the relationship between mass and metabolism, which can be described by the power function y = a*x^b
ggplot(ds1, aes(x = log_mass, y = SMR, group = ID_fish)) + # plot raw SMR data by ID against log10(mass)
geom_point(size = 2) + # size of data points
facet_wrap(. ~ Temp_class) + # produces one plot per temperature treatment
theme_bw() +
scale_y_log10() # plot data with log10 scaled y axis
The trends become linear when SMR is presented on log scale
Univariate SMR models for two temperature treatments
# model structure using log10 SMR and mass for fish at 18C only
SMR18_mod<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish), # random slopes by month, random intercepts by ID_fish
data=ds18, REML = T)
summary(SMR18_mod) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
## Data: ds18
##
## REML criterion at convergence: -553.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6860 -0.7401 0.0462 0.7148 2.7433
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 1.063e-17 3.261e-09
## Residual 9.832e-03 9.916e-02
## Number of obs: 317, groups: ID_fish, 72
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.67878 0.01415 315.00000 -47.96 <2e-16 ***
## log_mass 0.78058 0.02562 315.00000 30.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log_mass -0.919
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')plot(SMR18_mod) # plot model residuals
ranova(SMR18_mod) # anova-like test shows insignificant random slopes## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 276.78 -545.56
## (1 | ID_fish) 3 276.78 -547.56 0 1 1SMR18_mod_1<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish), # random slopes term removed
data=ds18, REML = T)
summary(SMR18_mod_1) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
## Data: ds18
##
## REML criterion at convergence: -553.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6860 -0.7401 0.0462 0.7148 2.7433
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 1.063e-17 3.261e-09
## Residual 9.832e-03 9.916e-02
## Number of obs: 317, groups: ID_fish, 72
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.67878 0.01415 315.00000 -47.96 <2e-16 ***
## log_mass 0.78058 0.02562 315.00000 30.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log_mass -0.919
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')plot(SMR18_mod_1) # plot model residuals
ranova(SMR18_mod_1) # anova-like test shows insignificant random intercepts, but retain in model for testing## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 276.78 -545.56
## (1 | ID_fish) 3 276.78 -547.56 0 1 1# model structure using log10 SMR and mass for fish at 23C only
SMR23_mod<-lmer(log10(SMR) ~ log_mass + (1 + Month_int|ID_fish), # random slopes by month, random intercepts by ID_fish
data=ds23, REML=T)
summary(SMR23_mod) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 + Month_int | ID_fish)
## Data: ds23
##
## REML criterion at convergence: -649.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1396 -0.5503 -0.0023 0.5686 3.2950
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.0052305 0.07232
## Month_int 0.0003171 0.01781 -0.94
## Residual 0.0056041 0.07486
## Number of obs: 307, groups: ID_fish, 64
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.47671 0.01133 101.03860 -42.07 <2e-16 ***
## log_mass 0.61389 0.02319 128.47109 26.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log_mass -0.876plot(SMR23_mod) # plot model residuals
ranova(SMR23_mod) # anova-like test shows insignificant random slopes## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + (1 + Month_int | ID_fish)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 6 324.70 -637.41
## Month_int in (1 + Month_int | ID_fish) 4 320.96 -633.92 7.4863 2 0.02368
##
## <none>
## Month_int in (1 + Month_int | ID_fish) *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1SMR23_mod_1<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish), # random slopes term removed
data=ds23, REML=T)
summary(SMR23_mod_1) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
## Data: ds23
##
## REML criterion at convergence: -641.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0810 -0.5308 -0.0445 0.5962 3.1721
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.0006479 0.02545
## Residual 0.0064034 0.08002
## Number of obs: 307, groups: ID_fish, 64
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.47203 0.01064 133.59553 -44.38 <2e-16 ***
## log_mass 0.60664 0.02244 178.85157 27.03 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log_mass -0.851plot(SMR23_mod_1) # plot model residuals
ranova(SMR23_mod_1) # anova-like test shows insignificant random intercepts, but retain in model for testing## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 320.96 -633.92
## (1 | ID_fish) 3 318.78 -631.56 4.3593 1 0.03681 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Combined SMR model
# model using log10 SMR and the interaction of log10 mass and temperature treatment
SMR_mod <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Sex*as.factor(Temp_class) + Sex*as.factor(Temp_class)*log_mass + Density + (1|ID_fish),
data=ds1, REML=T)
plot(SMR_mod) # generate model output
summary(SMR_mod) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Sex * as.factor(Temp_class) +
## Sex * as.factor(Temp_class) * log_mass + Density + (1 | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -832.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.87616 -0.67366 0.01528 0.66654 2.90845
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.0001859 0.01364
## Residual 0.0082037 0.09057
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.676339 0.031184 156.099873 -21.689
## log_mass 0.796149 0.032692 255.489558 24.353
## as.factor(Temp_class)23 0.165816 0.033347 188.402172 4.972
## SexM -0.057479 0.032974 278.565445 -1.743
## Density -0.002566 0.003298 152.321130 -0.778
## log_mass:as.factor(Temp_class)23 -0.089839 0.063547 212.044147 -1.414
## as.factor(Temp_class)23:SexM 0.086503 0.047344 220.930300 1.827
## log_mass:SexM 0.160075 0.063848 346.693411 2.507
## log_mass:as.factor(Temp_class)23:SexM -0.200908 0.093156 271.605997 -2.157
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 1.48e-06 ***
## SexM 0.0824 .
## Density 0.4377
## log_mass:as.factor(Temp_class)23 0.1589
## as.factor(Temp_class)23:SexM 0.0690 .
## log_mass:SexM 0.0126 *
## log_mass:as.factor(Temp_class)23:SexM 0.0319 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss as.(T_)23 SexM Densty lg_:.(T_)23 a.(T_)23: lg_:SM
## log_mass -0.663
## as.fc(T_)23 -0.476 0.565
## SexM -0.372 0.559 0.359
## Density -0.769 0.091 0.122 -0.018
## lg_:.(T_)23 0.390 -0.520 -0.925 -0.286 -0.111
## a.(T_)23:SM 0.294 -0.393 -0.698 -0.696 -0.032 0.646
## log_mss:SxM 0.290 -0.506 -0.282 -0.925 0.018 0.259 0.644
## l_:.(T_)23: -0.231 0.351 0.626 0.633 0.030 -0.677 -0.916 -0.685ranova(SMR_mod) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class) + as.factor(Temp_class):Sex + log_mass:Sex + log_mass:as.factor(Temp_class):Sex
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 416.12 -810.24
## (1 | ID_fish) 10 415.93 -811.85 0.38676 1 0.534# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
SMR_mod_month <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),
data=ds1, REML=T)
plot(SMR_mod_month) # generate model output
summary(SMR_mod_month) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 + Month_int | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -839.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.95365 -0.62500 0.02573 0.65621 2.95222
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 1.150e-03 0.033906
## Month_int 4.241e-05 0.006512 -1.00
## Residual 8.165e-03 0.090361
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.693448 0.029255 209.328902 -23.704
## log_mass 0.835361 0.027904 309.422980 29.937
## as.factor(Temp_class)23 0.195188 0.023171 174.966869 8.424
## SexM 0.014368 0.009421 121.002750 1.525
## Density -0.003163 0.003287 220.153668 -0.962
## log_mass:as.factor(Temp_class)23 -0.154606 0.043270 289.652921 -3.573
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 1.29e-14 ***
## SexM 0.129854
## Density 0.336940
## log_mass:as.factor(Temp_class)23 0.000413 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexM Densty
## log_mass -0.619
## as.fc(T_)23 -0.481 0.639
## SexM -0.212 0.215 0.100
## Density -0.816 0.113 0.125 -0.020
## lg_:.(T_)23 0.432 -0.651 -0.915 -0.145 -0.112
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')ranova(SMR_mod_month) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 10 419.78 -819.56
## Month_int in (1 + Month_int | ID_fish) 8 419.39 -822.78 0.78507 2
## Pr(>Chisq)
## <none>
## Month_int in (1 + Month_int | ID_fish) 0.6753# model with random slopes by temp treatment
SMR_mod_trt <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0+Temp_class|ID_fish),
data=ds1, REML=T)
plot(SMR_mod_trt) # generate model output
summary(SMR_mod_trt) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (0 + Temp_class | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -838.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.95727 -0.63762 0.00871 0.65226 2.96919
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish Temp_class18 4.911e-05 0.007008
## Temp_class23 1.882e-04 0.013718 0.56
## Residual 8.348e-03 0.091368
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.696017 0.028900 153.536942 -24.084
## log_mass 0.834470 0.027427 187.550606 30.425
## as.factor(Temp_class)23 0.195939 0.022670 184.796780 8.643
## SexM 0.015418 0.009284 87.658721 1.661
## Density -0.002715 0.003269 150.116946 -0.831
## log_mass:as.factor(Temp_class)23 -0.154162 0.042967 211.007690 -3.588
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 2.56e-15 ***
## SexM 0.100349
## Density 0.407416
## log_mass:as.factor(Temp_class)23 0.000414 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexM Densty
## log_mass -0.609
## as.fc(T_)23 -0.471 0.632
## SexM -0.211 0.218 0.096
## Density -0.823 0.112 0.126 -0.020
## lg_:.(T_)23 0.420 -0.643 -0.912 -0.141 -0.108ranova(SMR_mod_trt) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 10 419.42 -818.84
## Temp_class in (0 + Temp_class | ID_fish) 8 419.39 -822.78 0.061499 2
## Pr(>Chisq)
## <none>
## Temp_class in (0 + Temp_class | ID_fish) 0.9697# model with random intercepts by temp treatment
SMR_mod_trtint <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),
data=ds1, REML=T)
plot(SMR_mod_trtint) # generate model output
summary(SMR_mod_trtint) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 | Temp_class)
## Data: ds1
##
## REML criterion at convergence: -838.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.96092 -0.63658 0.00346 0.65601 2.98231
##
## Random effects:
## Groups Name Variance Std.Dev.
## Temp_class (Intercept) 0.0002785 0.01669
## Residual 0.0084475 0.09191
## Number of obs: 453, groups: Temp_class, 2
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -6.951e-01 3.302e-02 2.625e-08 -21.051
## log_mass 8.339e-01 2.736e-02 4.470e+02 30.485
## as.factor(Temp_class)23 1.941e-01 3.247e-02 6.139e-09 5.978
## SexM 1.524e-02 9.057e-03 4.470e+02 1.682
## Density -2.806e-03 3.211e-03 4.470e+02 -0.874
## log_mass:as.factor(Temp_class)23 -1.502e-01 4.233e-02 4.470e+02 -3.548
## Pr(>|t|)
## (Intercept) 1.000000
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 1.000000
## SexM 0.093181 .
## Density 0.382747
## log_mass:as.factor(Temp_class)23 0.000429 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexM Densty
## log_mass -0.528
## as.fc(T_)23 -0.545 0.440
## SexM -0.180 0.216 0.067
## Density -0.707 0.109 0.088 -0.021
## lg_:.(T_)23 0.371 -0.651 -0.629 -0.141 -0.111ranova(SMR_mod_trtint) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 419.33 -822.66
## (1 | Temp_class) 7 419.33 -824.66 -1.1369e-13 1 1# test sex terms for model inclusion
SMR_mod_sex <- lmer(log10(SMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass + Density + (1 |ID_fish),
data=ds1, REML=T)
plot(SMR_mod_sex)
summary(SMR_mod_sex)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log10(SMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass +
## Density + (1 | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -832.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.87616 -0.67366 0.01528 0.66654 2.90845
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.0001859 0.01364
## Residual 0.0082037 0.09057
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error
## Density -0.002566 0.003298
## as.factor(Sex)F:as.factor(Temp_class)18 -0.676339 0.031184
## as.factor(Sex)M:as.factor(Temp_class)18 -0.733818 0.035973
## as.factor(Sex)F:as.factor(Temp_class)23 -0.510523 0.033096
## as.factor(Sex)M:as.factor(Temp_class)23 -0.481499 0.030654
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass 0.796149 0.032692
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass 0.956224 0.055066
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass 0.706309 0.054270
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass 0.665476 0.040820
## df t value Pr(>|t|)
## Density 152.321130 -0.778 0.438
## as.factor(Sex)F:as.factor(Temp_class)18 156.099871 -21.689 <2e-16
## as.factor(Sex)M:as.factor(Temp_class)18 226.413501 -20.399 <2e-16
## as.factor(Sex)F:as.factor(Temp_class)23 149.558833 -15.426 <2e-16
## as.factor(Sex)M:as.factor(Temp_class)23 158.947247 -15.707 <2e-16
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass 255.489556 24.353 <2e-16
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass 374.137705 17.365 <2e-16
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass 198.971185 13.015 <2e-16
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass 256.884341 16.303 <2e-16
##
## Density
## as.factor(Sex)F:as.factor(Temp_class)18 ***
## as.factor(Sex)M:as.factor(Temp_class)18 ***
## as.factor(Sex)F:as.factor(Temp_class)23 ***
## as.factor(Sex)M:as.factor(Temp_class)23 ***
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## Densty as.(S)F:.(T_)18 as.(S)M:.(T_)18 as.(S)F:.(T_)23
## as.(S)F:.(T_)18 -0.769
## as.(S)M:.(T_)18 -0.683 0.526
## as.(S)F:.(T_)23 -0.602 0.463 0.411
## as.(S)M:.(T_)23 -0.719 0.553 0.491 0.433
## a.(S)F:.(T_)18: 0.091 -0.663 -0.062 -0.055
## a.(S)M:.(T_)18: 0.074 -0.057 -0.729 -0.045
## a.(S)F:.(T_)23: -0.075 0.058 0.051 -0.694
## a.(S)M:.(T_)23: -0.004 0.003 0.003 0.002
## as.(S)M:.(T_)23 a.(S)F:.(T_)18: a.(S)M:.(T_)18: a.(S)F:.(T_)23:
## as.(S)F:.(T_)18
## as.(S)M:.(T_)18
## as.(S)F:.(T_)23
## as.(S)M:.(T_)23
## a.(S)F:.(T_)18: -0.065
## a.(S)M:.(T_)18: -0.054 0.007
## a.(S)F:.(T_)23: 0.054 -0.007 -0.006
## a.(S)M:.(T_)23: -0.604 0.000 0.000 0.000Again, our ranova analysis shows that the random intercepts by ID_fish are not significant to this model. However, this will be retained because repeated measures were made on the same individuals, and so that the model can be included in the below covariance matrix.
Plot SMR vs. log mass mean predicted level
pred_SMR <- ggpredict(SMR_mod, terms = c("log_mass", "Temp_class"))
ggplot() +
geom_errorbar(data = pred_SMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) + # add error bars showing 95% confidence intervals
geom_point(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group)) +
geom_line(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group)) +
labs(x = "log_mass", y = "Predicted SMR", color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10() 
F test on interaction effect
anova(SMR_mod, type = "3")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## log_mass 9.2265 9.2265 1 269.79 1124.6757
## as.factor(Temp_class) 0.6283 0.6283 1 222.17 76.5912
## Sex 0.0030 0.0030 1 223.64 0.3604
## Density 0.0050 0.0050 1 152.32 0.6056
## log_mass:as.factor(Temp_class) 0.1351 0.1351 1 273.42 16.4634
## as.factor(Temp_class):Sex 0.0274 0.0274 1 220.93 3.3384
## log_mass:Sex 0.0134 0.0134 1 274.33 1.6352
## log_mass:as.factor(Temp_class):Sex 0.0382 0.0382 1 271.61 4.6513
## Pr(>F)
## log_mass < 2.2e-16 ***
## as.factor(Temp_class) 5.351e-16 ***
## Sex 0.54888
## Density 0.43767
## log_mass:as.factor(Temp_class) 6.477e-05 ***
## as.factor(Temp_class):Sex 0.06903 .
## log_mass:Sex 0.20207
## log_mass:as.factor(Temp_class):Sex 0.03191 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Generate fitted data
z_SMR<-augment(SMR_mod) # get variables and fitted values from model
head(z_SMR)## # A tibble: 6 × 18
## .rownames `log10(SMR)` log_mass `as.factor(Temp_class)` Sex Density ID_fish
## <chr> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
## 1 1 -0.0733 0.423 23 F 8 81
## 2 2 -0.546 -0.0915 23 M 8 82
## 3 3 -0.276 0.250 23 F 8 83
## 4 4 -0.0948 0.436 23 M 8 84
## 5 5 -0.151 0.360 23 M 8 85
## 6 7 -0.296 0.310 23 M 8 87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## # .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## # .weights <dbl>, .wtres <dbl>z_SMR <- z_SMR %>%
mutate(
pwr.fitted = 10^(.fitted)
)
##### check assumptions: plot all residual values against the fitted
ggplot(z_SMR, aes(x = .fitted, y = .resid)) + # same plot as before
geom_point(size = 2) 
ggplot(z_SMR, aes(x = log_mass, y = SMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`), # note this will look odd if there are additional fixed effects
alpha = 0.6) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10()+
labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))
At this stage, we have found a strong model for SMR which captures variance by treatment and uses centered mass. The next step is to repeat this modeling approach for RMR and MMR.
RMR MODEL
# model using log10 RMR and the interaction of log10 mass and temperature treatment
RMR_mod <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Sex*as.factor(Temp_class) + Sex*as.factor(Temp_class)*log_mass + Density + (1 |ID_fish),
data=ds1, REML=T)
plot(RMR_mod)
summary(RMR_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Sex * as.factor(Temp_class) +
## Sex * as.factor(Temp_class) * log_mass + Density + (1 | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -895.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6217 -0.6035 0.0362 0.6597 2.4226
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.001370 0.03701
## Residual 0.006309 0.07943
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.458863 0.033569 211.189136 -13.669
## log_mass 0.715302 0.032984 359.999351 21.686
## as.factor(Temp_class)23 0.133553 0.035094 249.693441 3.806
## SexM -0.053158 0.033240 333.370196 -1.599
## Density -0.004299 0.003555 213.356425 -1.209
## log_mass:as.factor(Temp_class)23 -0.097402 0.065508 312.941087 -1.487
## as.factor(Temp_class)23:SexM 0.070119 0.049047 273.943840 1.430
## log_mass:SexM 0.124388 0.062014 428.861974 2.006
## log_mass:as.factor(Temp_class)23:SexM -0.163282 0.093374 373.092786 -1.749
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 0.000178 ***
## SexM 0.110717
## Density 0.227872
## log_mass:as.factor(Temp_class)23 0.138054
## as.factor(Temp_class)23:SexM 0.153960
## log_mass:SexM 0.045505 *
## log_mass:as.factor(Temp_class)23:SexM 0.081167 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss as.(T_)23 SexM Densty lg_:.(T_)23 a.(T_)23: lg_:SM
## log_mass -0.655
## as.fc(T_)23 -0.456 0.542
## SexM -0.381 0.555 0.362
## Density -0.791 0.134 0.124 0.004
## lg_:.(T_)23 0.363 -0.509 -0.904 -0.280 -0.110
## a.(T_)23:SM 0.289 -0.382 -0.710 -0.678 -0.042 0.642
## log_mss:SxM 0.295 -0.523 -0.280 -0.895 -0.003 0.263 0.607
## l_:.(T_)23: -0.224 0.352 0.629 0.595 0.038 -0.697 -0.887 -0.664ranova(RMR_mod)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class) + as.factor(Temp_class):Sex + log_mass:Sex + log_mass:as.factor(Temp_class):Sex
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 447.95 -873.91
## (1 | ID_fish) 10 437.96 -855.92 19.988 1 7.792e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
RMR_mod_month <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),
data=ds1, REML=T)
plot(RMR_mod_month) # generate model output
summary(RMR_mod_month) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 + Month_int | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -904.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6681 -0.5983 0.0284 0.6278 2.4239
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 1.509e-03 0.038848
## Month_int 1.135e-06 0.001065 -1.00
## Residual 6.380e-03 0.079873
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.476881 0.031184 208.946377 -15.293
## log_mass 0.748619 0.027720 371.916742 27.006
## as.factor(Temp_class)23 0.159088 0.023714 198.905581 6.709
## SexM 0.003107 0.011023 86.886610 0.282
## Density -0.004258 0.003515 203.921547 -1.211
## log_mass:as.factor(Temp_class)23 -0.153474 0.043242 358.751139 -3.549
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 2e-10 ***
## SexM 0.778752
## Density 0.227163
## log_mass:as.factor(Temp_class)23 0.000438 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexM Densty
## log_mass -0.610
## as.fc(T_)23 -0.458 0.615
## SexM -0.203 0.178 0.064
## Density -0.836 0.152 0.128 -0.012
## lg_:.(T_)23 0.404 -0.644 -0.886 -0.118 -0.112
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00514734 (tol = 0.002, component 1)ranova(RMR_mod_month) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 10 452.28 -884.56
## Month_int in (1 + Month_int | ID_fish) 8 452.24 -888.49 0.076141 2
## Pr(>Chisq)
## <none>
## Month_int in (1 + Month_int | ID_fish) 0.9626# model with random intercepts by temp treatment
RMR_mod_trtint <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),
data=ds1, REML=T)
plot(RMR_mod_trtint)
summary(RMR_mod_trtint)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 | Temp_class)
## Data: ds1
##
## REML criterion at convergence: -886.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0205 -0.6163 0.0171 0.6597 3.1772
##
## Random effects:
## Groups Name Variance Std.Dev.
## Temp_class (Intercept) 9.358e-05 0.009674
## Residual 7.588e-03 0.087111
## Number of obs: 453, groups: Temp_class, 2
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -4.775e-01 2.868e-02 3.908e-08 -16.649
## log_mass 7.495e-01 2.593e-02 4.470e+02 28.908
## as.factor(Temp_class)23 1.463e-01 2.518e-02 5.802e-09 5.811
## SexM 3.679e-03 8.584e-03 4.470e+02 0.429
## Density -4.302e-03 3.044e-03 4.470e+02 -1.413
## log_mass:as.factor(Temp_class)23 -1.232e-01 4.012e-02 4.470e+02 -3.071
## Pr(>|t|)
## (Intercept) 1.00000
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 1.00000
## SexM 0.66846
## Density 0.15824
## log_mass:as.factor(Temp_class)23 0.00226 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexM Densty
## log_mass -0.576
## as.fc(T_)23 -0.510 0.538
## SexM -0.197 0.216 0.081
## Density -0.772 0.109 0.108 -0.021
## lg_:.(T_)23 0.404 -0.651 -0.768 -0.141 -0.111ranova(RMR_mod_trtint)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 443.3 -870.6
## (1 | Temp_class) 7 443.3 -872.6 2.2737e-13 1 1# model with random slopes by temp treatment
RMR_mod_trt <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0 + Temp_class|ID_fish),
data=ds1, REML=T)
plot(RMR_mod_trt)
summary(RMR_mod_trt)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (0 + Temp_class | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -904.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6323 -0.5883 0.0271 0.6401 2.3369
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish Temp_class18 0.001392 0.03731
## Temp_class23 0.001059 0.03254 0.15
## Residual 0.006388 0.07992
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.476861 0.031190 208.937568 -15.289
## log_mass 0.748448 0.027931 302.207977 26.796
## as.factor(Temp_class)23 0.158264 0.023388 232.444603 6.767
## SexM 0.003190 0.010952 82.507699 0.291
## Density -0.004258 0.003507 196.713602 -1.214
## log_mass:as.factor(Temp_class)23 -0.150036 0.042951 298.007754 -3.493
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 1.06e-10 ***
## SexM 0.77155
## Density 0.22619
## log_mass:as.factor(Temp_class)23 0.00055 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexM Densty
## log_mass -0.607
## as.fc(T_)23 -0.470 0.624
## SexM -0.197 0.169 0.057
## Density -0.835 0.151 0.135 -0.013
## lg_:.(T_)23 0.410 -0.653 -0.886 -0.112 -0.117
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluesranova(RMR_mod_trt)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 10 452.33 -884.66
## Temp_class in (0 + Temp_class | ID_fish) 8 452.24 -888.49 0.17085 2
## Pr(>Chisq)
## <none>
## Temp_class in (0 + Temp_class | ID_fish) 0.9181# test sex terms for model inclusion
RMR_mod_sex <- lmer(log10(RMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass + Density + (1 |ID_fish),
data=ds1, REML=T)
plot(RMR_mod_sex)
summary(RMR_mod_sex)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log10(RMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass +
## Density + (1 | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -895.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6217 -0.6035 0.0362 0.6597 2.4226
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.001370 0.03701
## Residual 0.006309 0.07943
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error
## Density -0.004299 0.003555
## as.factor(Sex)F:as.factor(Temp_class)18 -0.458863 0.033569
## as.factor(Sex)M:as.factor(Temp_class)18 -0.512022 0.037156
## as.factor(Sex)F:as.factor(Temp_class)23 -0.325311 0.035831
## as.factor(Sex)M:as.factor(Temp_class)23 -0.308350 0.032984
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass 0.715302 0.032984
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass 0.839690 0.052870
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass 0.617900 0.056382
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass 0.579006 0.041101
## df t value Pr(>|t|)
## Density 213.356426 -1.209 0.228
## as.factor(Sex)F:as.factor(Temp_class)18 211.189136 -13.669 <2e-16
## as.factor(Sex)M:as.factor(Temp_class)18 283.322955 -13.780 <2e-16
## as.factor(Sex)F:as.factor(Temp_class)23 200.632088 -9.079 <2e-16
## as.factor(Sex)M:as.factor(Temp_class)23 207.520955 -9.348 <2e-16
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass 359.999351 21.686 <2e-16
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass 439.883035 15.882 <2e-16
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass 294.895487 10.959 <2e-16
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass 367.788434 14.088 <2e-16
##
## Density
## as.factor(Sex)F:as.factor(Temp_class)18 ***
## as.factor(Sex)M:as.factor(Temp_class)18 ***
## as.factor(Sex)F:as.factor(Temp_class)23 ***
## as.factor(Sex)M:as.factor(Temp_class)23 ***
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## Densty as.(S)F:.(T_)18 as.(S)M:.(T_)18 as.(S)F:.(T_)23
## as.(S)F:.(T_)18 -0.791
## as.(S)M:.(T_)18 -0.711 0.562
## as.(S)F:.(T_)23 -0.620 0.490 0.441
## as.(S)M:.(T_)23 -0.731 0.578 0.520 0.453
## a.(S)F:.(T_)18: 0.134 -0.655 -0.095 -0.083
## a.(S)M:.(T_)18: 0.080 -0.063 -0.687 -0.050
## a.(S)F:.(T_)23: -0.049 0.039 0.035 -0.681
## a.(S)M:.(T_)23: 0.015 -0.012 -0.011 -0.010
## as.(S)M:.(T_)23 a.(S)F:.(T_)18: a.(S)M:.(T_)18: a.(S)F:.(T_)23:
## as.(S)F:.(T_)18
## as.(S)M:.(T_)18
## as.(S)F:.(T_)23
## as.(S)M:.(T_)23
## a.(S)F:.(T_)18: -0.098
## a.(S)M:.(T_)18: -0.059 0.011
## a.(S)F:.(T_)23: 0.036 -0.007 -0.004
## a.(S)M:.(T_)23: -0.579 0.002 0.001 -0.001# plot RMR vs. log mass mean predicted level
pred_RMR <- ggpredict(RMR_mod, terms = c("log_mass", "Temp_class"))
ggplot() +
geom_errorbar(data = pred_RMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) + # add error bars showing 95% confidence intervals
geom_point(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group)) +
geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group)) +
labs(x = "log_mass_centre", y = "Predicted RMR", color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10() 
# plot the raw data as a cloud around the predictions
ggplot() +
geom_point(data = ds1, aes(x = log_mass, y = RMR, group = factor(Temp_class), color = factor(Temp_class)),
alpha = 0.3, position = position_jitter(width = 0.1)) +
geom_point(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group),
size = 3) +
geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group),
size = 2) +
labs(x = "log_mass_centre", y = "Predicted RMR", color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10() 
# F test on interaction effect
anova(RMR_mod, type = "3")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## log_mass 5.4611 5.4611 1 367.11 865.6183
## as.factor(Temp_class) 0.2934 0.2934 1 276.53 46.5051
## Sex 0.0034 0.0034 1 276.43 0.5449
## Density 0.0092 0.0092 1 213.36 1.4625
## log_mass:as.factor(Temp_class) 0.0917 0.0917 1 375.49 14.5343
## as.factor(Temp_class):Sex 0.0129 0.0129 1 273.94 2.0439
## log_mass:Sex 0.0053 0.0053 1 375.33 0.8386
## log_mass:as.factor(Temp_class):Sex 0.0193 0.0193 1 373.09 3.0579
## Pr(>F)
## log_mass < 2.2e-16 ***
## as.factor(Temp_class) 5.738e-11 ***
## Sex 0.4610311
## Density 0.2278721
## log_mass:as.factor(Temp_class) 0.0001608 ***
## as.factor(Temp_class):Sex 0.1539602
## log_mass:Sex 0.3603912
## log_mass:as.factor(Temp_class):Sex 0.0811666 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# generate fitted data
z_RMR<-augment(RMR_mod) # get variables and fitted values from model
head(z_RMR)## # A tibble: 6 × 18
## .rownames `log10(RMR)` log_mass `as.factor(Temp_class)` Sex Density ID_fish
## <chr> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
## 1 1 -0.0212 0.423 23 F 8 81
## 2 2 -0.386 -0.0915 23 M 8 82
## 3 3 -0.147 0.250 23 F 8 83
## 4 4 0.00674 0.436 23 M 8 84
## 5 5 -0.0573 0.360 23 M 8 85
## 6 7 -0.161 0.310 23 M 8 87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## # .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## # .weights <dbl>, .wtres <dbl>z_RMR <- z_RMR %>%
mutate(
pwr.fitted = 10^(.fitted)
)
# check assumptions: plot all residual values against the fitted
ggplot(z_RMR, aes(x = .fitted, y = .resid)) + # same plot as before
geom_point(size = 2) 
ggplot(z_RMR, aes(x = log_mass, y = RMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),
alpha = 0.6) +
theme_bw() +
scale_y_log10() +
scale_color_manual(values = c("#78B7C5", "#F21A00"))
MMR MODEL
# model using log10 MMR and the interaction of log10 mass and temperature treatment
MMR_mod <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 |ID_fish),
data=ds1, REML=T)
plot(MMR_mod)
summary(MMR_mod) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -1086.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0358 -0.6335 0.0066 0.5706 3.5127
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.0008517 0.02918
## Residual 0.0042430 0.06514
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.034542 0.025478 210.005598 -1.356
## log_mass 0.718537 0.022624 372.713216 31.760
## as.factor(Temp_class)23 0.050788 0.019226 264.403760 2.642
## SexM 0.026480 0.009007 86.996070 2.940
## Density -0.012575 0.002877 211.708210 -4.371
## log_mass:as.factor(Temp_class)23 -0.096076 0.035291 359.553907 -2.722
## Pr(>|t|)
## (Intercept) 0.17663
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 0.00874 **
## SexM 0.00421 **
## Density 1.94e-05 ***
## log_mass:as.factor(Temp_class)23 0.00680 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexM Densty
## log_mass -0.605
## as.fc(T_)23 -0.457 0.614
## SexM -0.202 0.176 0.060
## Density -0.838 0.152 0.130 -0.013
## lg_:.(T_)23 0.400 -0.644 -0.884 -0.115 -0.111ranova(MMR_mod)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 543.31 -1070.6
## (1 | ID_fish) 7 534.20 -1054.4 18.232 1 1.956e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
MMR_mod_month <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),
data=ds1, REML=T)
plot(MMR_mod_month) # generate model output
summary(MMR_mod_month) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 + Month_int | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -1090.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6553 -0.6206 0.0243 0.5554 3.5381
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.0013142 0.03625
## Month_int 0.0001224 0.01106 -0.62
## Residual 0.0039368 0.06274
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.037512 0.025652 199.690054 -1.462
## log_mass 0.716517 0.023329 271.438607 30.714
## as.factor(Temp_class)23 0.051810 0.019205 169.950435 2.698
## SexM 0.028424 0.008949 87.097271 3.176
## Density -0.012195 0.002938 191.627985 -4.151
## log_mass:as.factor(Temp_class)23 -0.096354 0.036605 238.653992 -2.632
## Pr(>|t|)
## (Intercept) 0.14523
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 0.00768 **
## SexM 0.00206 **
## Density 4.99e-05 ***
## log_mass:as.factor(Temp_class)23 0.00903 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexM Densty
## log_mass -0.576
## as.fc(T_)23 -0.449 0.613
## SexM -0.195 0.178 0.055
## Density -0.840 0.116 0.123 -0.019
## lg_:.(T_)23 0.382 -0.639 -0.885 -0.108 -0.093ranova(MMR_mod_month) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 10 545.27 -1070.5
## Month_int in (1 + Month_int | ID_fish) 8 543.31 -1070.6 3.9032 2 0.142# model with random intercepts by temp treatment
MMR_mod_trtint <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),
data=ds1, REML=T)
plot(MMR_mod_trtint)
summary(MMR_mod_trtint) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 | Temp_class)
## Data: ds1
##
## REML criterion at convergence: -1068.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3156 -0.6143 0.0195 0.6269 3.1035
##
## Random effects:
## Groups Name Variance Std.Dev.
## Temp_class (Intercept) 6.744e-05 0.008212
## Residual 5.053e-03 0.071082
## Number of obs: 453, groups: Temp_class, 2
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -5.223e-02 2.351e-02 1.327e-07 -2.221
## log_mass 7.324e-01 2.116e-02 4.470e+02 34.618
## as.factor(Temp_class)23 5.136e-02 2.079e-02 2.029e-08 2.470
## SexM 2.738e-02 7.005e-03 4.470e+02 3.909
## Density -1.111e-02 2.484e-03 4.470e+02 -4.473
## log_mass:as.factor(Temp_class)23 -9.468e-02 3.274e-02 4.470e+02 -2.892
## Pr(>|t|)
## (Intercept) 0.999999
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 1.000000
## SexM 0.000107 ***
## Density 9.79e-06 ***
## log_mass:as.factor(Temp_class)23 0.004014 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexM Densty
## log_mass -0.574
## as.fc(T_)23 -0.512 0.532
## SexM -0.196 0.216 0.080
## Density -0.768 0.109 0.107 -0.021
## lg_:.(T_)23 0.402 -0.651 -0.759 -0.141 -0.111ranova(MMR_mod_trtint)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 534.2 -1052.4
## (1 | Temp_class) 7 534.2 -1054.4 6.8212e-13 1 1# model with random slopes by temp treatment
MMR_mod_trt <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0 + Temp_class|ID_fish),
data=ds1, REML=T)
plot(MMR_mod_trt)
summary(MMR_mod_trt) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (0 + Temp_class | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -1086.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0147 -0.6251 0.0075 0.5777 3.5188
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish Temp_class18 0.0008129 0.02851
## Temp_class23 0.0009152 0.03025 0.43
## Residual 0.0042424 0.06513
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.034292 0.025461 208.686591 -1.347
## log_mass 0.718885 0.022540 289.137042 31.894
## as.factor(Temp_class)23 0.051319 0.019279 235.372611 2.662
## SexM 0.026620 0.009019 86.512834 2.952
## Density -0.012647 0.002880 209.059535 -4.391
## log_mass:as.factor(Temp_class)23 -0.097241 0.035385 315.337174 -2.748
## Pr(>|t|)
## (Intercept) 0.17950
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 0.00831 **
## SexM 0.00407 **
## Density 1.79e-05 ***
## log_mass:as.factor(Temp_class)23 0.00634 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexM Densty
## log_mass -0.604
## as.fc(T_)23 -0.451 0.609
## SexM -0.204 0.180 0.062
## Density -0.839 0.151 0.128 -0.012
## lg_:.(T_)23 0.396 -0.639 -0.883 -0.117 -0.109ranova(MMR_mod_trt)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 10 543.33 -1066.7
## Temp_class in (0 + Temp_class | ID_fish) 8 543.31 -1070.6 0.035729 2
## Pr(>Chisq)
## <none>
## Temp_class in (0 + Temp_class | ID_fish) 0.9823# test sex terms for model inclusion
MMR_mod_sex <- lmer(log10(MMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass + Density + (1 |ID_fish),
data=ds1, REML=T)
plot(MMR_mod_sex)
summary(MMR_mod_sex)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log10(MMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass +
## Density + (1 | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -1073.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1037 -0.6119 0.0047 0.5638 3.5177
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.0008985 0.02997
## Residual 0.0042347 0.06507
## Number of obs: 453, groups: ID_fish, 91
##
## Fixed effects:
## Estimate Std. Error
## Density -0.012730 0.002902
## as.factor(Sex)F:as.factor(Temp_class)18 -0.041773 0.027406
## as.factor(Sex)M:as.factor(Temp_class)18 0.008437 0.030355
## as.factor(Sex)F:as.factor(Temp_class)23 0.006331 0.029249
## as.factor(Sex)M:as.factor(Temp_class)23 0.048888 0.026927
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass 0.731205 0.026963
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass 0.688049 0.043257
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass 0.650875 0.046066
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass 0.606438 0.033600
## df t value Pr(>|t|)
## Density 211.516807 -4.386 1.82e-05
## as.factor(Sex)F:as.factor(Temp_class)18 209.572960 -1.524 0.1290
## as.factor(Sex)M:as.factor(Temp_class)18 281.936327 0.278 0.7812
## as.factor(Sex)F:as.factor(Temp_class)23 199.134068 0.216 0.8289
## as.factor(Sex)M:as.factor(Temp_class)23 206.145308 1.816 0.0709
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass 357.870749 27.119 < 2e-16
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass 439.391714 15.906 < 2e-16
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass 292.475337 14.129 < 2e-16
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass 365.600428 18.049 < 2e-16
##
## Density ***
## as.factor(Sex)F:as.factor(Temp_class)18
## as.factor(Sex)M:as.factor(Temp_class)18
## as.factor(Sex)F:as.factor(Temp_class)23
## as.factor(Sex)M:as.factor(Temp_class)23 .
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## Densty as.(S)F:.(T_)18 as.(S)M:.(T_)18 as.(S)F:.(T_)23
## as.(S)F:.(T_)18 -0.790
## as.(S)M:.(T_)18 -0.711 0.562
## as.(S)F:.(T_)23 -0.620 0.490 0.440
## as.(S)M:.(T_)23 -0.731 0.578 0.520 0.453
## a.(S)F:.(T_)18: 0.133 -0.656 -0.095 -0.083
## a.(S)M:.(T_)18: 0.080 -0.063 -0.688 -0.050
## a.(S)F:.(T_)23: -0.049 0.039 0.035 -0.682
## a.(S)M:.(T_)23: 0.015 -0.012 -0.011 -0.009
## as.(S)M:.(T_)23 a.(S)F:.(T_)18: a.(S)M:.(T_)18: a.(S)F:.(T_)23:
## as.(S)F:.(T_)18
## as.(S)M:.(T_)18
## as.(S)F:.(T_)23
## as.(S)M:.(T_)23
## a.(S)F:.(T_)18: -0.097
## a.(S)M:.(T_)18: -0.059 0.011
## a.(S)F:.(T_)23: 0.036 -0.007 -0.004
## a.(S)M:.(T_)23: -0.580 0.002 0.001 -0.001# plot MMR vs. log mass mean predicted level
pred_MMR <- ggpredict(MMR_mod, terms = c("log_mass", "Temp_class"))
ggplot() +
geom_errorbar(data = pred_MMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) + # add error bars showing 95% confidence intervals
geom_point(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group)) +
geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group)) +
labs(x = "log_mass", y = "Predicted MMR", color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10() 
# Plot the raw data as a cloud around the predictions
ggplot() +
geom_point(data = ds1, aes(x = log_mass, y = MMR, group = factor(Temp_class), color = factor(Temp_class)),
alpha = 0.3, position = position_jitter(width = 0.1)) +
geom_point(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group),
size = 3) +
geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group),
size = 2) +
labs(x = "log_mass_centre", y = "Predicted MMR", color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10() 
# F test on interaction effect
anova(MMR_mod, type = "3")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_mass 6.1674 6.1674 1 349.13 1453.5364 < 2.2e-16
## as.factor(Temp_class) 0.0296 0.0296 1 264.40 6.9784 0.008741
## Sex 0.0367 0.0367 1 87.00 8.6421 0.004206
## Density 0.0811 0.0811 1 211.71 19.1076 1.936e-05
## log_mass:as.factor(Temp_class) 0.0314 0.0314 1 359.55 7.4113 0.006797
##
## log_mass ***
## as.factor(Temp_class) **
## Sex **
## Density ***
## log_mass:as.factor(Temp_class) **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Generate fitted data
z_MMR<-augment(MMR_mod) # get variables and fitted values from model
head(z_MMR)## # A tibble: 6 × 18
## .rownames `log10(MMR)` log_mass `as.factor(Temp_class)` Sex Density ID_fish
## <chr> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
## 1 1 0.154 0.423 23 F 8 81
## 2 2 0.0458 -0.0915 23 M 8 82
## 3 3 0.159 0.250 23 F 8 83
## 4 4 0.215 0.436 23 M 8 84
## 5 5 0.227 0.360 23 M 8 85
## 6 7 0.0919 0.310 23 M 8 87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## # .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## # .weights <dbl>, .wtres <dbl>z_MMR <- z_MMR %>%
mutate(
pwr.fitted = 10^(.fitted)
)
# check assumptions: plot all residual values against the fitted
ggplot(z_MMR, aes(x = .fitted, y = .resid)) + # same plot as before
geom_point(size = 2) 
ggplot(z_MMR, aes(x = log_mass, y = MMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),
alpha = 0.6) +
theme_bw() +
scale_y_log10()+
scale_color_manual(values = c("#78B7C5", "#F21A00"))
FIGURE 2
Generate a figure that displays mass-scaling relationships for all three metabolic rates at both temperatures
### FIG. 2 ###
fig2.1 <- ggplot() +
geom_ribbon(data = pred_SMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) + # add error bars showing 95% confidence intervals
geom_line(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group),
size = 1) +
geom_point(data = ds1, aes(x = log_mass, y = SMR, group = factor(Temp_class), color = factor(Temp_class)),
alpha = 0.25, position = position_jitter(width = 0.1)) +
labs(x = "log10(mass (g))", y = expression(SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
xlim(-0.25, 1.1) +
scale_y_log10(limits = c(0.15, 1.6))
fig2.2 <- ggplot() +
geom_ribbon(data = pred_RMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) + # add error bars showing 95% confidence intervals
geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group),
size = 1) +
geom_point(data = ds1, aes(x = log_mass, y = RMR, group = factor(Temp_class), color = factor(Temp_class)),
alpha = 0.25, position = position_jitter(width = 0.1)) +
labs(x = "log10(mass (g))", y = expression(RMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
xlim(-0.25, 1.1) +
scale_y_log10(limits=c(0.25,2.1))
fig2.3 <- ggplot() +
geom_ribbon(data = pred_MMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) + # add error bars showing 95% confidence intervals
geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group),
size = 1) +
geom_point(data = ds1, aes(x = log_mass, y = MMR, group = factor(Temp_class), color = factor(Temp_class)),
alpha = 0.25, position = position_jitter(width = 0.1)) +
labs(x = expression(log[10]~(mass~(g))), y = expression(MMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
xlim(-0.25, 1.1) +
scale_y_log10(limits=c(0.5, 5))
fig2.4 <- ggplot(z_SMR, aes(x = log_mass, y = SMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`), # note this will look odd if there are additional fixed effects
alpha = 0.6) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10(limits = c(0.15, 1.6))+
xlim(-0.25, 1.1) +
labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))
fig2.5 <- ggplot(z_RMR, aes(x = log_mass, y = RMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`), # note this will look odd if there are additional fixed effects
alpha = 0.6) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10(limits=c(0.25,2.1))+
xlim(-0.25, 1.1) +
labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~RMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))
fig2.6 <- ggplot(z_MMR, aes(x = log_mass, y = MMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`), # note this will look odd if there are additional fixed effects
alpha = 0.6) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10(limits=c(0.5, 5))+
xlim(-0.25, 1.1) +
labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~MMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))
library(ggpubr)
fig2 <- ggarrange(print(fig2.1 + rremove("xlab")),
print(fig2.2 + rremove("xlab")),
print(fig2.3 + rremove("xlab")),
print(fig2.4 + rremove("xlab") + rremove("ylab")),
print(fig2.5 + rremove("xlab") + rremove("ylab")),
print(fig2.6 + rremove("xlab") + rremove("ylab")),
ncol = 3, nrow = 2, widths = c(5, 5, 5),
align = "v",
common.legend = TRUE,
legend = "right")
require(grid)
annotate_figure(fig2, bottom = textGrob(expression(log[10]~(mass~(g))), gp = gpar(cex = 1.1)))
setwd(figures_wd)
jpeg(file = "Fig2.jpeg", width = 1440, height=960, units = "px")
annotate_figure(fig2, bottom = textGrob(expression(log[10]~(mass~(g))), gp = gpar(cex = 1.1)))
dev.off()## png
## 2UNIVARIATE MODEL OF FAT MASS
#subset ds1 to only include fish for which fat mass was measured
ds_fat <- ds1[which(!is.na(ds1$log_fatmass) & ds1$Month == "September"),]
fat_mod <- lm(log_fatmass ~ Sex + Sex:log_finalmass_centre + Sex:Temp_class + Temp_class + log_finalmass_centre:Temp_class + Sex:Temp_class:log_finalmass_centre + Density, data=ds_fat)
fat_mod1 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + SMR:Temp_class + Sex + Density, data=ds_fat) #using only final measurement of SMR here
fat_mod2 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + RMR:Temp_class + Sex + Density, data=ds_fat) #using only final measurement of RMR here
fat_mod3 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + MMR:Temp_class + Sex + Density, data=ds_fat) #using only final measurement of MMR here
summary(fat_mod) # generate model output##
## Call:
## lm(formula = log_fatmass ~ Sex + Sex:log_finalmass_centre + Sex:Temp_class +
## Temp_class + log_finalmass_centre:Temp_class + Sex:Temp_class:log_finalmass_centre +
## Density, data = ds_fat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.66014 -0.10148 0.03335 0.10829 0.39864
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.633537 0.145266 -11.245 < 2e-16
## SexM -0.284274 0.179915 -1.580 0.118
## Temp_class23 -0.207258 0.198303 -1.045 0.299
## Density -0.008147 0.013416 -0.607 0.545
## SexF:log_finalmass_centre 1.480710 0.186207 7.952 9.34e-12
## SexM:log_finalmass_centre 2.108583 0.246715 8.547 6.26e-13
## SexM:Temp_class23 -0.124996 0.272443 -0.459 0.648
## log_finalmass_centre:Temp_class23 0.616740 0.357173 1.727 0.088
## SexM:log_finalmass_centre:Temp_class23 -0.030108 0.491759 -0.061 0.951
##
## (Intercept) ***
## SexM
## Temp_class23
## Density
## SexF:log_finalmass_centre ***
## SexM:log_finalmass_centre ***
## SexM:Temp_class23
## log_finalmass_centre:Temp_class23 .
## SexM:log_finalmass_centre:Temp_class23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1789 on 81 degrees of freedom
## (23 observations deleted due to missingness)
## Multiple R-squared: 0.807, Adjusted R-squared: 0.7879
## F-statistic: 42.34 on 8 and 81 DF, p-value: < 2.2e-16plot(fat_mod) # residuals plot
anova(fat_mod)## Analysis of Variance Table
##
## Response: log_fatmass
## Df Sum Sq Mean Sq F value Pr(>F)
## Sex 1 0.2138 0.2138 6.6797 0.01155 *
## Temp_class 1 0.0241 0.0241 0.7533 0.38801
## Density 1 0.0041 0.0041 0.1276 0.72191
## Sex:log_finalmass_centre 2 10.3111 5.1556 161.0605 < 2e-16 ***
## Sex:Temp_class 1 0.0995 0.0995 3.1072 0.08172 .
## log_finalmass_centre:Temp_class 1 0.1894 0.1894 5.9183 0.01719 *
## Sex:log_finalmass_centre:Temp_class 1 0.0001 0.0001 0.0037 0.95133
## Residuals 81 2.5928 0.0320
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Marginal Effects
#Plotting marginal effects using the code Pete gave me
require(ggeffects)
pred.fat_mod <- ggpredict(fat_mod, terms = c("log_finalmass_centre", "Sex", "Temp_class"))
plot(pred.fat_mod)
ggplot() +
geom_line(data = pred.fat_mod, aes(x = x, y = predicted, color = group, group = group)) +
geom_ribbon(data = pred.fat_mod, aes(x= x, ymin = conf.low, ymax = conf.high, fill = group, group = group), alpha = 0.4) +
labs(x = expression(log[10]~(centred~mass~(g))), y = expression(log[10]~(fat~mass~(g))), color = expression("Sex")) +
theme_classic() +
scale_color_manual(values = c("darkgreen", "purple")) +
scale_fill_manual(values = c("darkgreen", "purple")) 
summary(fat_mod1) # generate model output##
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class +
## SMR:Temp_class + Sex + Density, data = ds_fat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75613 -0.09036 0.04041 0.12172 0.31448
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.669290 0.132302 -12.617 < 2e-16 ***
## Temp_class23 -0.335509 0.143164 -2.344 0.0215 *
## SexM -0.028106 0.045131 -0.623 0.5352
## Density -0.009403 0.014587 -0.645 0.5210
## Temp_class18:log_finalmass_centre 1.373244 0.306634 4.478 2.41e-05 ***
## Temp_class23:log_finalmass_centre 2.827861 0.397690 7.111 3.91e-10 ***
## Temp_class18:SMR 0.136718 0.146222 0.935 0.3525
## Temp_class23:SMR -0.268150 0.236012 -1.136 0.2592
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1866 on 82 degrees of freedom
## (23 observations deleted due to missingness)
## Multiple R-squared: 0.7876, Adjusted R-squared: 0.7694
## F-statistic: 43.43 on 7 and 82 DF, p-value: < 2.2e-16plot(fat_mod) # residuals plot
summary(fat_mod2) # generate model output##
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class +
## RMR:Temp_class + Sex + Density, data = ds_fat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.76919 -0.08987 0.04334 0.12239 0.30783
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.67164 0.13345 -12.526 < 2e-16 ***
## Temp_class23 -0.33243 0.14187 -2.343 0.0215 *
## SexM -0.02392 0.04425 -0.540 0.5904
## Density -0.01111 0.01425 -0.780 0.4379
## Temp_class18:log_finalmass_centre 1.40379 0.31899 4.401 3.22e-05 ***
## Temp_class23:log_finalmass_centre 2.74098 0.41355 6.628 3.35e-09 ***
## Temp_class18:RMR 0.09706 0.12531 0.775 0.4409
## Temp_class23:RMR -0.15521 0.18365 -0.845 0.4005
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1876 on 82 degrees of freedom
## (23 observations deleted due to missingness)
## Multiple R-squared: 0.7853, Adjusted R-squared: 0.7669
## F-statistic: 42.84 on 7 and 82 DF, p-value: < 2.2e-16plot(fat_mod) # residuals plot
summary(fat_mod3) # generate model output##
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class +
## MMR:Temp_class + Sex + Density, data = ds_fat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75804 -0.09095 0.02879 0.12219 0.33877
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.666135 0.133622 -12.469 < 2e-16 ***
## Temp_class23 -0.320143 0.141473 -2.263 0.0263 *
## SexM -0.009131 0.042876 -0.213 0.8319
## Density -0.014221 0.014104 -1.008 0.3163
## Temp_class18:log_finalmass_centre 1.787541 0.296595 6.027 4.57e-08 ***
## Temp_class23:log_finalmass_centre 2.718573 0.301131 9.028 6.33e-14 ***
## Temp_class18:MMR -0.035603 0.058075 -0.613 0.5415
## Temp_class23:MMR -0.075744 0.061506 -1.231 0.2217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1869 on 82 degrees of freedom
## (23 observations deleted due to missingness)
## Multiple R-squared: 0.7867, Adjusted R-squared: 0.7685
## F-statistic: 43.21 on 7 and 82 DF, p-value: < 2.2e-16plot(fat_mod) # residuals plot
pred_fat <- ggpredict(fat_mod, terms = c("log_finalmass_centre", "Temp_class"))
# Plot the raw data
ggplot() +
geom_point(data = ds_fat, aes(x = log_finalmass_centre, y = log_fatmass, group = factor(Temp_class), color = factor(Temp_class))) +
geom_line(data = pred_fat, aes(x = x, y = predicted, color = group, group = group)) +
geom_ribbon(data = pred_fat, aes(x= x, ymin = conf.low, ymax = conf.high, fill = group, group = group), alpha = 0.4) +
labs(x = expression(log[10]~(centred~mass~(g))), y = expression(log[10]~(fat~mass~(g))), color = expression("Temperature ("*degree*C*")")) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_fill_manual(values = c("#78B7C5", "#F21A00")) 
ggplot(data = ds_fat, aes(y = Fat_mass_g, x = Sex, fill = factor(Sex))) +
scale_fill_manual(values = c("darkgreen", "purple")) +
geom_boxplot() +
labs(x = "Sex", y = "Fat mass (g)") +
theme_classic()
# F test on interaction term
anova(fat_mod)## Analysis of Variance Table
##
## Response: log_fatmass
## Df Sum Sq Mean Sq F value Pr(>F)
## Sex 1 0.2138 0.2138 6.6797 0.01155 *
## Temp_class 1 0.0241 0.0241 0.7533 0.38801
## Density 1 0.0041 0.0041 0.1276 0.72191
## Sex:log_finalmass_centre 2 10.3111 5.1556 161.0605 < 2e-16 ***
## Sex:Temp_class 1 0.0995 0.0995 3.1072 0.08172 .
## log_finalmass_centre:Temp_class 1 0.1894 0.1894 5.9183 0.01719 *
## Sex:log_finalmass_centre:Temp_class 1 0.0001 0.0001 0.0037 0.95133
## Residuals 81 2.5928 0.0320
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fat_mod1)## Analysis of Variance Table
##
## Response: log_fatmass
## Df Sum Sq Mean Sq F value Pr(>F)
## Temp_class 1 0.0406 0.0406 1.1650 0.28359
## Sex 1 0.1974 0.1974 5.6706 0.01957 *
## Density 1 0.0041 0.0041 0.1173 0.73285
## Temp_class:log_finalmass_centre 2 10.2619 5.1309 147.4086 < 2e-16 ***
## Temp_class:SMR 2 0.0769 0.0384 1.1045 0.33625
## Residuals 82 2.8542 0.0348
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fat_mod2)## Analysis of Variance Table
##
## Response: log_fatmass
## Df Sum Sq Mean Sq F value Pr(>F)
## Temp_class 1 0.0406 0.0406 1.1526 0.28616
## Sex 1 0.1974 0.1974 5.6100 0.02021 *
## Density 1 0.0041 0.0041 0.1161 0.73423
## Temp_class:log_finalmass_centre 2 10.2619 5.1309 145.8320 < 2e-16 ***
## Temp_class:RMR 2 0.0460 0.0230 0.6542 0.52257
## Residuals 82 2.8851 0.0352
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fat_mod3)## Analysis of Variance Table
##
## Response: log_fatmass
## Df Sum Sq Mean Sq F value Pr(>F)
## Temp_class 1 0.0406 0.0406 1.1604 0.28455
## Sex 1 0.1974 0.1974 5.6480 0.01981 *
## Density 1 0.0041 0.0041 0.1168 0.73337
## Temp_class:log_finalmass_centre 2 10.2619 5.1309 146.8201 < 2e-16 ***
## Temp_class:MMR 2 0.0654 0.0327 0.9364 0.39619
## Residuals 82 2.8657 0.0349
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This model is excluded from the following covariance matrix because we were only able to collect one measurement of fat mass from each fish, all at the very end of the experiment. Without an initial measurement of fat mass (logistically impossible as it was a lethal measurement), we are unable to assess within-subject residual correlation. For this reason, we have assessed it separately.
COVARIANCE MATRIX
covariance matrix with only sex specific traits
RN <- bf(mass1 ~ 0 + Sex:Temp_class + Sex:Temp_class:Month_centred_factor + Density + (1 + Month_centred |a| ID_fish)) +
bf(SMR1 ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 |a| ID_fish)) +
bf(RMR1 ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 |a| ID_fish)) +
bf(MMR1 ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 |a| ID_fish)) +
set_rescor(T)
prior1 <- c(set_prior("normal(0,2)", class = "b", resp = c("mass1", "SMR1", "RMR1","MMR1")),
set_prior("cauchy(0,2)", class = "sd", resp = c("mass1", "SMR1", "RMR1","MMR1")),
set_prior("lkj(2)", class = "cor"))
fit0 <- brm(formula = RN,
data = ds1,
prior = prior1,
seed = 42, warmup = 500, iter = 6000, chains = 4, cores = 4,
control=list(adapt_delta=0.97, max_treedepth = 15),
save_ranef = T
)
print(fit0, digits = 3)## Family: MV(gaussian, gaussian, gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: mass1 ~ 0 + Sex:Temp_class + Sex:Temp_class:Month_centred_factor + Density + (1 + Month_centred | a | ID_fish)
## SMR1 ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 | a | ID_fish)
## RMR1 ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 | a | ID_fish)
## MMR1 ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 | a | ID_fish)
## Data: ds1 (Number of observations: 453)
## Draws: 4 chains, each with iter = 6000; warmup = 500; thin = 1;
## total post-warmup draws = 22000
##
## Multilevel Hyperparameters:
## ~ID_fish (Number of levels: 91)
## Estimate Est.Error l-95% CI u-95% CI
## sd(mass1_Intercept) 0.167 0.014 0.142 0.197
## sd(mass1_Month_centred) 0.039 0.004 0.032 0.047
## sd(SMR1_Intercept) 0.012 0.007 0.001 0.026
## sd(RMR1_Intercept) 0.032 0.005 0.023 0.043
## sd(MMR1_Intercept) 0.031 0.005 0.021 0.041
## cor(mass1_Intercept,mass1_Month_centred) -0.358 0.103 -0.545 -0.145
## cor(mass1_Intercept,SMR1_Intercept) -0.110 0.304 -0.667 0.512
## cor(mass1_Month_centred,SMR1_Intercept) 0.331 0.286 -0.341 0.787
## cor(mass1_Intercept,RMR1_Intercept) 0.160 0.152 -0.149 0.444
## cor(mass1_Month_centred,RMR1_Intercept) 0.179 0.151 -0.121 0.464
## cor(SMR1_Intercept,RMR1_Intercept) 0.394 0.316 -0.375 0.832
## cor(mass1_Intercept,MMR1_Intercept) 0.179 0.167 -0.161 0.488
## cor(mass1_Month_centred,MMR1_Intercept) 0.199 0.153 -0.111 0.488
## cor(SMR1_Intercept,MMR1_Intercept) 0.008 0.308 -0.595 0.598
## cor(RMR1_Intercept,MMR1_Intercept) 0.399 0.165 0.051 0.694
## Rhat Bulk_ESS Tail_ESS
## sd(mass1_Intercept) 1.000 5430 9939
## sd(mass1_Month_centred) 1.000 5852 9991
## sd(SMR1_Intercept) 1.001 3809 8610
## sd(RMR1_Intercept) 1.001 2154 6458
## sd(MMR1_Intercept) 1.000 7495 11259
## cor(mass1_Intercept,mass1_Month_centred) 1.000 6044 10623
## cor(mass1_Intercept,SMR1_Intercept) 1.000 15888 15169
## cor(mass1_Month_centred,SMR1_Intercept) 1.000 17212 14454
## cor(mass1_Intercept,RMR1_Intercept) 1.000 11506 14354
## cor(mass1_Month_centred,RMR1_Intercept) 1.001 10437 15802
## cor(SMR1_Intercept,RMR1_Intercept) 1.003 1116 2049
## cor(mass1_Intercept,MMR1_Intercept) 1.000 12958 15452
## cor(mass1_Month_centred,MMR1_Intercept) 1.000 13864 16602
## cor(SMR1_Intercept,MMR1_Intercept) 1.001 2043 4458
## cor(RMR1_Intercept,MMR1_Intercept) 1.000 7519 14396
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI
## mass1_Density 0.012 0.006 0.000
## mass1_SexF:Temp_class18 0.284 0.054 0.180
## mass1_SexM:Temp_class18 0.206 0.060 0.088
## mass1_SexF:Temp_class23 0.230 0.057 0.119
## mass1_SexM:Temp_class23 0.148 0.060 0.031
## mass1_SexF:Temp_class18:Month_centred_factor1 0.161 0.016 0.130
## mass1_SexM:Temp_class18:Month_centred_factor1 0.157 0.021 0.117
## mass1_SexF:Temp_class23:Month_centred_factor1 0.104 0.021 0.062
## mass1_SexM:Temp_class23:Month_centred_factor1 0.161 0.023 0.116
## mass1_SexF:Temp_class18:Month_centred_factor2 0.254 0.020 0.214
## mass1_SexM:Temp_class18:Month_centred_factor2 0.174 0.027 0.121
## mass1_SexF:Temp_class23:Month_centred_factor2 0.164 0.026 0.114
## mass1_SexM:Temp_class23:Month_centred_factor2 0.270 0.028 0.215
## mass1_SexF:Temp_class18:Month_centred_factor3 0.198 0.026 0.147
## mass1_SexM:Temp_class18:Month_centred_factor3 0.114 0.034 0.047
## mass1_SexF:Temp_class23:Month_centred_factor3 0.169 0.033 0.105
## mass1_SexM:Temp_class23:Month_centred_factor3 0.277 0.036 0.207
## mass1_SexF:Temp_class18:Month_centred_factor4 0.407 0.031 0.346
## mass1_SexM:Temp_class18:Month_centred_factor4 0.340 0.040 0.260
## mass1_SexF:Temp_class23:Month_centred_factor4 0.288 0.040 0.209
## mass1_SexM:Temp_class23:Month_centred_factor4 0.425 0.044 0.339
## SMR1_Density -0.002 0.003 -0.009
## SMR1_SexF:Temp_class18 -0.875 0.040 -0.953
## SMR1_SexM:Temp_class18 -0.967 0.048 -1.062
## SMR1_SexF:Temp_class23 -0.683 0.045 -0.772
## SMR1_SexM:Temp_class23 -0.647 0.040 -0.726
## SMR1_SexF:Temp_class18:log_mass_centre 0.806 0.037 0.734
## SMR1_SexM:Temp_class18:log_mass_centre 0.963 0.058 0.850
## SMR1_SexF:Temp_class23:log_mass_centre 0.710 0.058 0.599
## SMR1_SexM:Temp_class23:log_mass_centre 0.674 0.045 0.584
## RMR1_Density -0.006 0.003 -0.012
## RMR1_SexF:Temp_class18 -0.615 0.039 -0.691
## RMR1_SexM:Temp_class18 -0.689 0.045 -0.779
## RMR1_SexF:Temp_class23 -0.466 0.045 -0.554
## RMR1_SexM:Temp_class23 -0.436 0.039 -0.513
## RMR1_SexF:Temp_class18:log_mass_centre 0.710 0.035 0.642
## RMR1_SexM:Temp_class18:log_mass_centre 0.821 0.053 0.718
## RMR1_SexF:Temp_class23:log_mass_centre 0.622 0.056 0.514
## RMR1_SexM:Temp_class23:log_mass_centre 0.579 0.042 0.494
## MMR1_Density -0.013 0.003 -0.019
## MMR1_SexF:Temp_class18 -0.204 0.035 -0.272
## MMR1_SexM:Temp_class18 -0.148 0.040 -0.225
## MMR1_SexF:Temp_class23 -0.132 0.040 -0.210
## MMR1_SexM:Temp_class23 -0.085 0.035 -0.153
## MMR1_SexF:Temp_class18:log_mass_centre 0.718 0.031 0.657
## MMR1_SexM:Temp_class18:log_mass_centre 0.678 0.045 0.590
## MMR1_SexF:Temp_class23:log_mass_centre 0.629 0.049 0.533
## MMR1_SexM:Temp_class23:log_mass_centre 0.593 0.037 0.519
## u-95% CI Rhat Bulk_ESS Tail_ESS
## mass1_Density 0.024 1.000 6711 10981
## mass1_SexF:Temp_class18 0.390 1.000 4564 8643
## mass1_SexM:Temp_class18 0.324 1.000 4386 8415
## mass1_SexF:Temp_class23 0.343 1.000 4265 7980
## mass1_SexM:Temp_class23 0.268 1.000 4726 7882
## mass1_SexF:Temp_class18:Month_centred_factor1 0.193 1.001 11148 14768
## mass1_SexM:Temp_class18:Month_centred_factor1 0.198 1.000 10838 15271
## mass1_SexF:Temp_class23:Month_centred_factor1 0.145 1.000 11252 14898
## mass1_SexM:Temp_class23:Month_centred_factor1 0.205 1.000 10625 14688
## mass1_SexF:Temp_class18:Month_centred_factor2 0.294 1.001 5942 11371
## mass1_SexM:Temp_class18:Month_centred_factor2 0.227 1.002 5705 10164
## mass1_SexF:Temp_class23:Month_centred_factor2 0.216 1.000 6792 11809
## mass1_SexM:Temp_class23:Month_centred_factor2 0.326 1.000 6631 9866
## mass1_SexF:Temp_class18:Month_centred_factor3 0.248 1.001 4768 8353
## mass1_SexM:Temp_class18:Month_centred_factor3 0.181 1.002 4628 9215
## mass1_SexF:Temp_class23:Month_centred_factor3 0.235 1.000 5556 9519
## mass1_SexM:Temp_class23:Month_centred_factor3 0.348 1.000 5297 8763
## mass1_SexF:Temp_class18:Month_centred_factor4 0.468 1.001 4489 8492
## mass1_SexM:Temp_class18:Month_centred_factor4 0.419 1.002 4337 8685
## mass1_SexF:Temp_class23:Month_centred_factor4 0.367 1.000 5037 9083
## mass1_SexM:Temp_class23:Month_centred_factor4 0.510 1.000 4961 7621
## SMR1_Density 0.004 1.000 12077 16800
## SMR1_SexF:Temp_class18 -0.797 1.001 7008 11864
## SMR1_SexM:Temp_class18 -0.874 1.001 7995 12588
## SMR1_SexF:Temp_class23 -0.595 1.001 8792 12493
## SMR1_SexM:Temp_class23 -0.570 1.001 7953 13161
## SMR1_SexF:Temp_class18:log_mass_centre 0.878 1.001 8564 13224
## SMR1_SexM:Temp_class18:log_mass_centre 1.076 1.000 9974 14180
## SMR1_SexF:Temp_class23:log_mass_centre 0.824 1.001 10179 14311
## SMR1_SexM:Temp_class23:log_mass_centre 0.760 1.000 9639 14921
## RMR1_Density 0.001 1.000 10974 15126
## RMR1_SexF:Temp_class18 -0.537 1.001 7075 12036
## RMR1_SexM:Temp_class18 -0.601 1.001 7864 12329
## RMR1_SexF:Temp_class23 -0.379 1.001 8311 11983
## RMR1_SexM:Temp_class23 -0.359 1.000 7899 13195
## RMR1_SexF:Temp_class18:log_mass_centre 0.778 1.001 8538 13846
## RMR1_SexM:Temp_class18:log_mass_centre 0.927 1.000 10273 13978
## RMR1_SexF:Temp_class23:log_mass_centre 0.733 1.000 9957 13423
## RMR1_SexM:Temp_class23:log_mass_centre 0.661 1.000 9832 14947
## MMR1_Density -0.007 1.000 14150 16707
## MMR1_SexF:Temp_class18 -0.135 1.000 9836 14603
## MMR1_SexM:Temp_class18 -0.070 1.000 11607 14923
## MMR1_SexF:Temp_class23 -0.054 1.000 11757 14948
## MMR1_SexM:Temp_class23 -0.016 1.000 11091 15576
## MMR1_SexF:Temp_class18:log_mass_centre 0.779 1.000 12482 15761
## MMR1_SexM:Temp_class18:log_mass_centre 0.765 1.000 16329 16198
## MMR1_SexF:Temp_class23:log_mass_centre 0.725 1.000 15410 15767
## MMR1_SexM:Temp_class23:log_mass_centre 0.665 1.000 16232 16397
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_mass1 0.060 0.003 0.055 0.065 1.000 13626 15029
## sigma_SMR1 0.091 0.003 0.085 0.097 1.000 15574 16610
## sigma_RMR1 0.080 0.003 0.075 0.087 1.000 12035 17014
## sigma_MMR1 0.066 0.003 0.061 0.071 1.000 18093 16365
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(mass1,SMR1) -0.086 0.079 -0.239 0.070 1.000 8393 13542
## rescor(mass1,RMR1) -0.108 0.074 -0.252 0.037 1.000 8805 13667
## rescor(SMR1,RMR1) 0.862 0.013 0.834 0.886 1.000 16325 17291
## rescor(mass1,MMR1) -0.050 0.069 -0.185 0.086 1.000 14591 16962
## rescor(SMR1,MMR1) 0.143 0.050 0.046 0.239 1.001 26459 16631
## rescor(RMR1,MMR1) 0.207 0.050 0.108 0.303 1.001 26753 17238
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(fit0)
conditional_effects(fit0)
build covariance matrix with brms
RN <- bf(mass1 ~ 0 + Temp_class + Month_centred_factor:Temp_class + Sex + Density + (1 + Month_centred |a| ID_fish)) +
bf(SMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) +
bf(RMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) +
bf(MMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) +
set_rescor(T)
prior1 <- c(set_prior("normal(0,2)", class = "b", resp = c("mass1", "SMR1", "RMR1","MMR1")),
set_prior("cauchy(0,2)", class = "sd", resp = c("mass1", "SMR1", "RMR1","MMR1")),
set_prior("lkj(2)", class = "cor"))
fit1 <- brm(formula = RN,
data = ds1,
prior = prior1,
seed = 42, warmup = 500, iter = 10000, chains = 4, cores = 4,
control=list(adapt_delta=0.97, max_treedepth = 15),
save_ranef = T
)
print(fit1, digits = 3)## Family: MV(gaussian, gaussian, gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: mass1 ~ 0 + Temp_class + Month_centred_factor:Temp_class + Sex + Density + (1 + Month_centred | a | ID_fish)
## SMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish)
## RMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish)
## MMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish)
## Data: ds1 (Number of observations: 453)
## Draws: 4 chains, each with iter = 10000; warmup = 500; thin = 1;
## total post-warmup draws = 38000
##
## Multilevel Hyperparameters:
## ~ID_fish (Number of levels: 91)
## Estimate Est.Error l-95% CI u-95% CI
## sd(mass1_Intercept) 0.165 0.014 0.141 0.194
## sd(mass1_Month_centred) 0.040 0.004 0.034 0.048
## sd(SMR1_Intercept) 0.011 0.007 0.001 0.025
## sd(RMR1_Intercept) 0.031 0.005 0.023 0.042
## sd(MMR1_Intercept) 0.030 0.005 0.020 0.040
## cor(mass1_Intercept,mass1_Month_centred) -0.318 0.106 -0.512 -0.101
## cor(mass1_Intercept,SMR1_Intercept) -0.128 0.307 -0.688 0.503
## cor(mass1_Month_centred,SMR1_Intercept) 0.338 0.293 -0.346 0.800
## cor(mass1_Intercept,RMR1_Intercept) 0.158 0.151 -0.148 0.440
## cor(mass1_Month_centred,RMR1_Intercept) 0.177 0.149 -0.117 0.466
## cor(SMR1_Intercept,RMR1_Intercept) 0.363 0.330 -0.428 0.825
## cor(mass1_Intercept,MMR1_Intercept) 0.191 0.169 -0.156 0.505
## cor(mass1_Month_centred,MMR1_Intercept) 0.169 0.153 -0.137 0.460
## cor(SMR1_Intercept,MMR1_Intercept) 0.017 0.311 -0.592 0.610
## cor(RMR1_Intercept,MMR1_Intercept) 0.411 0.165 0.063 0.705
## Rhat Bulk_ESS Tail_ESS
## sd(mass1_Intercept) 1.001 8724 15843
## sd(mass1_Month_centred) 1.000 11231 18796
## sd(SMR1_Intercept) 1.001 7027 11325
## sd(RMR1_Intercept) 1.001 3711 10259
## sd(MMR1_Intercept) 1.000 13870 16558
## cor(mass1_Intercept,mass1_Month_centred) 1.001 9327 18688
## cor(mass1_Intercept,SMR1_Intercept) 1.000 33502 27685
## cor(mass1_Month_centred,SMR1_Intercept) 1.000 28395 23325
## cor(mass1_Intercept,RMR1_Intercept) 1.000 21696 25906
## cor(mass1_Month_centred,RMR1_Intercept) 1.000 16783 23670
## cor(SMR1_Intercept,RMR1_Intercept) 1.001 1755 3518
## cor(mass1_Intercept,MMR1_Intercept) 1.000 23420 26235
## cor(mass1_Month_centred,MMR1_Intercept) 1.000 26510 28843
## cor(SMR1_Intercept,MMR1_Intercept) 1.000 3062 7410
## cor(RMR1_Intercept,MMR1_Intercept) 1.000 13833 24495
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI
## mass1_Temp_class18 0.290 0.052 0.188 0.392
## mass1_Temp_class23 0.232 0.052 0.130 0.336
## mass1_SexM -0.072 0.034 -0.139 -0.004
## mass1_Density 0.011 0.006 -0.001 0.023
## mass1_Temp_class18:Month_centred_factor1 0.158 0.013 0.133 0.184
## mass1_Temp_class23:Month_centred_factor1 0.129 0.016 0.098 0.161
## mass1_Temp_class18:Month_centred_factor2 0.228 0.017 0.195 0.260
## mass1_Temp_class23:Month_centred_factor2 0.216 0.020 0.177 0.255
## mass1_Temp_class18:Month_centred_factor3 0.170 0.021 0.128 0.212
## mass1_Temp_class23:Month_centred_factor3 0.224 0.025 0.174 0.273
## mass1_Temp_class18:Month_centred_factor4 0.380 0.026 0.330 0.430
## mass1_Temp_class23:Month_centred_factor4 0.354 0.031 0.293 0.415
## SMR1_Temp_class18 -0.897 0.036 -0.969 -0.826
## SMR1_Temp_class23 -0.664 0.034 -0.731 -0.597
## SMR1_SexM 0.015 0.009 -0.004 0.034
## SMR1_Density -0.002 0.003 -0.009 0.004
## SMR1_Temp_class18:log_mass_centre 0.836 0.032 0.772 0.899
## SMR1_Temp_class23:log_mass_centre 0.681 0.037 0.608 0.753
## RMR1_Temp_class18 -0.631 0.036 -0.701 -0.561
## RMR1_Temp_class23 -0.446 0.034 -0.512 -0.378
## RMR1_SexM 0.002 0.010 -0.018 0.023
## RMR1_Density -0.006 0.003 -0.012 0.001
## RMR1_Temp_class18:log_mass_centre 0.732 0.030 0.673 0.791
## RMR1_Temp_class23:log_mass_centre 0.592 0.035 0.523 0.662
## MMR1_Temp_class18 -0.195 0.033 -0.259 -0.131
## MMR1_Temp_class23 -0.121 0.031 -0.182 -0.060
## MMR1_SexM 0.026 0.009 0.008 0.044
## MMR1_Density -0.013 0.003 -0.019 -0.007
## MMR1_Temp_class18:log_mass_centre 0.708 0.027 0.654 0.762
## MMR1_Temp_class23:log_mass_centre 0.609 0.032 0.547 0.672
## Rhat Bulk_ESS Tail_ESS
## mass1_Temp_class18 1.000 9281 16767
## mass1_Temp_class23 1.000 8639 15827
## mass1_SexM 1.000 5903 10692
## mass1_Density 1.000 14045 22047
## mass1_Temp_class18:Month_centred_factor1 1.000 18640 26859
## mass1_Temp_class23:Month_centred_factor1 1.000 18861 25964
## mass1_Temp_class18:Month_centred_factor2 1.000 10533 19970
## mass1_Temp_class23:Month_centred_factor2 1.000 11106 20122
## mass1_Temp_class18:Month_centred_factor3 1.000 8770 15803
## mass1_Temp_class23:Month_centred_factor3 1.000 9191 16982
## mass1_Temp_class18:Month_centred_factor4 1.000 8096 15500
## mass1_Temp_class23:Month_centred_factor4 1.000 8171 14915
## SMR1_Temp_class18 1.000 16774 22576
## SMR1_Temp_class23 1.000 16149 23956
## SMR1_SexM 1.000 45285 29965
## SMR1_Density 1.000 32156 28146
## SMR1_Temp_class18:log_mass_centre 1.000 19146 25494
## SMR1_Temp_class23:log_mass_centre 1.000 19297 25087
## RMR1_Temp_class18 1.000 14900 21747
## RMR1_Temp_class23 1.000 14417 22655
## RMR1_SexM 1.000 26137 29528
## RMR1_Density 1.000 25322 28582
## RMR1_Temp_class18:log_mass_centre 1.000 18250 25032
## RMR1_Temp_class23:log_mass_centre 1.000 17444 23199
## MMR1_Temp_class18 1.000 19049 24584
## MMR1_Temp_class23 1.000 17889 24390
## MMR1_SexM 1.000 30258 30267
## MMR1_Density 1.000 29008 30210
## MMR1_Temp_class18:log_mass_centre 1.000 23182 26563
## MMR1_Temp_class23:log_mass_centre 1.000 24019 27582
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_mass1 0.060 0.003 0.055 0.066 1.000 22767 28921
## sigma_SMR1 0.091 0.003 0.085 0.098 1.000 24535 28184
## sigma_RMR1 0.081 0.003 0.075 0.087 1.000 18866 28313
## sigma_MMR1 0.066 0.002 0.061 0.071 1.000 30665 26391
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(mass1,SMR1) -0.017 0.076 -0.166 0.133 1.000 17064 24877
## rescor(mass1,RMR1) -0.056 0.073 -0.198 0.087 1.000 17052 25013
## rescor(SMR1,RMR1) 0.863 0.013 0.836 0.887 1.000 27400 30731
## rescor(mass1,MMR1) -0.069 0.069 -0.203 0.065 1.000 26677 30302
## rescor(SMR1,MMR1) 0.138 0.049 0.040 0.234 1.000 45101 29153
## rescor(RMR1,MMR1) 0.203 0.049 0.105 0.298 1.000 45295 29787
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(fit1)
conditional_effects(fit1)
REPEATABILITY ESTIMATES
Use posterior distribution to generate repeatability estimates
View(posterior::summarise_draws(as_draws_array(fit1))) # view all predicted values of random effects
# random effects are fitted as standard dev, not variance, so we need to square them first
# var names to insert below: SMR, RMR, MMR, growth
# SMR repeatability
Vint_SMR <- (as_draws_array(fit1, variable = "sd_ID_fish__SMR1_Intercept"))^2
Vresid_SMR <- (as_draws_array(fit1, variable = "sigma_SMR1"))^2
R_SMR <- Vint_SMR / (Vint_SMR + Vresid_SMR)
"SMR"; mean(R_SMR); quantile(R_SMR, probs = c(0.025, 0.975))## [1] "SMR"## [1] 0.01979299## 2.5% 97.5%
## 5.327465e-05 7.128585e-02# RMR repeatability
Vint_RMR <- (as_draws_array(fit1, variable = "sd_ID_fish__RMR1_Intercept"))^2
Vresid_RMR <- (as_draws_array(fit1, variable = "sigma_RMR1"))^2
R_RMR <- Vint_RMR / (Vint_RMR + Vresid_RMR)
"RMR"; mean(R_RMR); quantile(R_RMR, probs = c(0.025, 0.975))## [1] "RMR"## [1] 0.130803## 2.5% 97.5%
## 0.06931255 0.21816931# MMR repeatability
Vint_MMR <- (as_draws_array(fit1, variable = "sd_ID_fish__MMR1_Intercept"))^2
Vresid_MMR <- (as_draws_array(fit1, variable = "sigma_MMR1"))^2
R_MMR <- Vint_MMR / (Vint_MMR + Vresid_MMR)
"MMR"; mean(R_MMR); quantile(R_MMR, probs = c(0.025, 0.975))## [1] "MMR"## [1] 0.1705791## 2.5% 97.5%
## 0.07647187 0.27647865# Growth rate repeatability
Vint_growth <- (as_draws_array(fit1, variable = "sd_ID_fish__mass1_Month_centred"))^2
Vslope_growth <- (as_draws_array(fit1, variable = "sd_ID_fish__mass1_Intercept"))^2
corr_growth <- (as_draws_array(fit1, variable = "cor_ID_fish__mass1_Intercept__mass1_Month_centred"))
Vresid_growth <- (as_draws_array(fit1, variable = "sigma_mass1"))^2
COVis_growth <- corr_growth*sqrt(Vslope_growth)*sqrt(Vint_growth)
x<- 0 #(can be 0 - 4)
VAR0 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R0_growth <- VAR0/(VAR0 + Vresid_growth)
"Growth - May"; mean(R0_growth); quantile(R0_growth, probs = c(0.025, 0.975))## [1] "Growth - May"## [1] 0.3126487## 2.5% 97.5%
## 0.2276433 0.4086757x<- 1 #(can be 0 - 4)
VAR1 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R1_growth <- VAR1/(VAR1 + Vresid_growth)
"Growth - June"; mean(R1_growth); quantile(R1_growth, probs = c(0.025, 0.975))## [1] "Growth - June"## [1] 0.8700228## 2.5% 97.5%
## 0.8252644 0.9079288x<- 2 #(can be 0 - 4)
VAR2 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R2_growth <- VAR2/(VAR2 + Vresid_growth)
"Growth - July"; mean(R2_growth); quantile(R2_growth, probs = c(0.025, 0.975))## [1] "Growth - July"## [1] 0.9651094## 2.5% 97.5%
## 0.9511292 0.9762125x<- 3 #(can be 0 - 4)
VAR3 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R3_growth <- VAR3/(VAR3 + Vresid_growth)
"Growth - August"; mean(R3_growth); quantile(R3_growth, probs = c(0.025, 0.975))## [1] "Growth - August"## [1] 0.9844639## 2.5% 97.5%
## 0.9780419 0.9894882x<- 4 #(can be 0 - 4)
VAR4 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R4_growth <- VAR4/(VAR4 + Vresid_growth)
"Growth - September"; mean(R4_growth); quantile(R4_growth, probs = c(0.025, 0.975))## [1] "Growth - September"## [1] 0.9912901## 2.5% 97.5%
## 0.9876578 0.9941299